program diag c ****************** Program DIAG ****************** c This program is designed to read in a real c symmetric matrix (but as a square matrix on c disk), diagonalize it, and return all c eigenvalues and corresponding eigenvectors. c c Copyright (C) 1992 Jeff Nichols and Jack Simons c c ************************************************** implicit real (a-h,o-z), integer(i-n) include "limits.h" c c Set core needed to: c 1 square array (eigenvectors) of order mmax c 1 symmetric array (input) of order mmax c 1 vector (eigenvalues) of length mmax c 1 scratch vector of length mmax c total = mmax*mmax + mmax*(mmax+1)/2 + 2*mmax c dimension pcore(mmaxsq + mmaxtr + 2*mmax), ipntr(10) itotal = mmax*mmax + mmax*(mmax+1)/2 + 2*mmax + 10 call modinfo('diag ') call shwmem('diag ',itotal) call goon 10 continue write(*,*) write(*,*)' Input the order of the real symmetric ' write(*,*)' array you would like diagonalized. ' read(*,*)msize if(msize.gt.mmax)then write(*,1010)mmax goto 10 endif ipntr(1) = 1 ipntr(2) = ipntr(1) + (msize*(msize+1))/2 ipntr(3) = ipntr(2) + msize*msize ipntr(4) = ipntr(3) + msize ipntr(5) = ipntr(4) + msize call driver(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)), & pcore(ipntr(4)),msize) call qmexit stop 1010 format(' Maximum order possible is less than ',I3,/, & ' Please try again. ') end subroutine driver(asym,avects,avals,scr1,msize) implicit real (a-h,o-z), integer(i-n) logical yesno character*30 dfname, fname, dscrpt character*1 ans dimension asym((msize*(msize+1))/2), & avects(msize,msize), & avals(msize), & scr1(msize) 10 continue dscrpt = 'matrix to be diagonalized' ids = 25 dfname = 'diag.in' idf = 7 fname = ' ' ifn = 30 iflnum = 32 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call rdmatd(fname,ifn,dscrpt,ids,avects,msize,msize, & iflnum) c c Close input file. c close(unit=32) else call rdmatu(dscrpt,ids,avects,msize,msize) endif c c Write input array back to user. c call wtmatu(dscrpt,ids,avects,msize,msize) write(*,*)' Is this correct? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then c c Pack the square array into a 1-dimensional symmetric c array. c ij = 0 do 40 i = 1,msize do 30 j = 1,i ij = ij + 1 asym(ij) = avects(i,j) avects(i,j) = 0.0e0 avects(j,i) = 0.0e0 30 continue 40 continue c c Diagonalize the array. c call diaghe(asym,msize,avals,msize,avects,msize,scr1) c c Write the eigenvalues to user. c write(*,*)' The eigenvalues of this matrix are: ' call matout(avals,msize,1,msize) c c Write the eigenvectors to user. c write(*,*)' their corresponding eigenvectors are: ' call matout(avects,msize,msize,msize) else write(*,*)' Please try again. ' goto 10 endif dscrpt = 'eigenvectors' ids = 12 dfname = 'diag_vectors.out' idf = 16 fname = ' ' ifn = 30 iflnum = 33 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call wtmatd(fname,ifn,dscrpt,ids,avects,msize,msize, & iflnum) c c Close output file. c close(unit=33) endif dscrpt = 'eigenvalues' ids = 11 dfname = 'diag_values.out' idf = 15 fname = ' ' ifn = 30 iflnum = 34 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call wtmatd(fname,ifn,dscrpt,ids,avals,msize,1, & iflnum) c c Close output file. c close(unit=34) endif return end