# Covariance Matrices Fortran Program

```        program qmask_demo
! Max Tegmark & Angelica de Oliveira-Costa 1999
! This code shows how to load a noise covariance matrix file.
! Must be compiled with real*8.
! For instance, on Solaris, you'd type f77 -r8 qmap_demo.f
!
! The matrices are stored in Solaris fortran binary form.
! Since they are symmetric, only the lower triangle is stored
! - the subroutines triangular2square and square2triangular
! convert between this and a regular square matrix format.
integer nmax, n, i, j
parameter(nmax=6500)
real M(nmax,nmax), sdev1, sdev2, cov, corr
character*60 fname
print *,'Number of pixels? (say 2695)'
if (n.gt.nmax) pause 'DEATH ERROR: nmax too small'
print *,'Name of covariance matrix file? (say 1Ka12_Sigma.raw)'
call triangular2square(nmax,n,M,M)
666     print *,'1st pixel number?'
print *,'2nd pixel number?'
if ((i.lt.1).or.(j.lt.1).or.(i.gt.n).or.(j.gt.n)) pause 'Bad pixel number'
sdev1   = sqrt(M(i,i))
sdev2   = sqrt(M(j,j))
cov     = M(i,j)
corr    = cov/(sdev1*sdev2)
print *,'RMS noise in 1st pixel.............',sdev1
print *,'RMS noise in 2nd pixel.............',sdev2
print *,'Noise covariance ..................',cov
print *,'Dimensionless noise correlation....',corr
goto 666
end

subroutine square2triangular(np,n,A,M)
! Converts the n x n symmetric matrix M into my 1-dimensional form.
! Because of the FORTRAN number ordering in matrices,
! it works even if A and M coincide physically.
implicit none
integer np, n, i, j
real A(np,np), M(n*(n+1)/2)
do i=1,n
do j=1,i
M((i*(i-1))/2 + j) = A(j,i)
end do
end do
return
end

subroutine triangular2square(np,n,M,A)
! The inverse of square2triangular.
! Unpacks the triangular M into the n x n symmetric matrix A.
! Because of the FORTRAN number ordering in matrices,
! it works even if A and M coincide physically.
implicit none
integer np, n, i, j
real A(np,np), M(n*(n+1)/2)
do i=n,1,-1
do j=i,1,-1
A(j,i) = M((i*(i-1))/2 + j)
end do
end do
do i=1,n
do j=1,i-1
A(i,j) = A(j,i)
end do
end do
return
end

subroutine SaveRawVector(n,a,fname)
integer n
real a(n)
character*60 fname
open(2,file=fname,form='unformatted')
print *,'Saving raw vector in ',fname
write(2) a
close(2)
print *,'File size should be ',8*n+8,' bytes.'
return
end

integer n
real a(n)
character*60 fname
open(2,file=fname,form='unformatted',status='old')