program matrixcossin c Matrix cosine and sine c ************************************************************ c n: Order of the matrix c tdelta: time step (decrease tdelta for higher accuracy) c tend: time $T$ upto which the time stepping is carried out c A(i,j): Input matrix (a(i,j) in the program) c program finds cos(AT) and sin(AT). c Redirect output to file: c Example matrixcossin.e >> matrixcossin.out c ************************************************************ integer i,j,k,m,n,nsteps real*8 tdelta,tend,a(100,100),asq(100,100),ident(100,100) complex*16 stiff(100,100),zn(100,100),znp1(100,100),b(100,100) open(unit=5,file='matrixcossin.dat',status='old') read(5,*) n,tdelta,tend do 10 i=1,n read(5,*) (a(i,j),j=1,n) 10 continue nsteps=nint(tend/tdelta) ident=0.0d0 asq=0.0d0 do 20 i=1,n ident(i,i)=1.0d0 zn(i,i)=1.0d0 20 continue do 40 i=1,n do 30 j=1,n asq(i,j)=0.0d0 do 25 k=1,n asq(i,j)=asq(i,j)+a(i,k)*a(k,j) 25 continue 30 continue 40 continue stiff(1:n,1:n)=6.0d0*(2.0d0*ident(1:n,1:n)-tdelta*(0.0d0,1.0d0)* & a(1:n,1:n))-tdelta**2*asq(1:n,1:n) b(1:n,1:n)=6.0d0*(2.0d0*ident(1:n,1:n)+tdelta*(0.0d0,1.0d0)* & a(1:n,1:n))-tdelta**2*asq(1:n,1:n) call zgaussj(stiff,n,100,b,n,100) do 100 m=1,nsteps do 70 i=1,n do 60 j=1,n znp1(i,j)=(0.0d0,0.0d0) do 50 k=1,n znp1(i,j)=znp1(i,j)+b(i,k)*zn(k,j) 50 continue 60 continue 70 continue zn(1:n,1:n)=znp1(1:n,1:n) 100 continue print*, 'Matrix cosine ' do 200 i=1,n print*, (dreal(znp1(i,j)),j=1,n) 200 continue print*,' ' print*, 'Matrix sine ' do 210 i=1,n print*, (dimag(znp1(i,j)),j=1,n) 210 continue end C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+ C File: zgaussj.f (Fortran 77 source for the zgaussj routine) C Author: Fredrik Jonsson C Date: December 26, 2005 C Last change: December 29, 2006 C Description: The zgaussj(a,n,np,b,m,mp) routine solves systems of C linear equations by Gauss-Jordan elimination, in similar C to as in Eq. (2.1.1) of W.H. Press et al, "Numerical C Recipes in Fortran 77, The Art of Scientific Computing", C 2nd Edn. (Cambridge University Press, Cambridge, 1992), C but here modified so as to solve complex-valued systems. C Here a(1:n,1:n) is the input system matrix stored in C an array of dimensions np by np, while b(1:n,1:m) is an C input matrix containing the m right-hand side vectors, C stored in an array of dimensions np by mp. On output, C a(1:n,1:n) is replaced by its matrix inverse, and C b(1:n,1:m) is replaced by the corresponding set of solu- C tion vectors. The parameter nmax contains the largest C anticipated value of n. C C Copyright (C) 2006, Fredrik Jonsson C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+ subroutine zgaussj(a,n,np,b,m,mp) integer m,mp,n,np,nmax complex*16 a(np,np),b(np,mp) parameter (nmax=400) integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax) complex*16 big,dum,pivinv irow=0 icol=0 do 11 j=1,n ipiv(j)=0 11 continue do 22 i=1,n big=0. do 13 j=1,n if(ipiv(j).ne.1)then do 12 k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.abs(big))then big=abs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then print*, 'Singular matrix ecountered in zgaussj' endif 12 continue endif 13 continue ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do 14 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 14 continue do 15 l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum 15 continue endif indxr(i)=irow indxc(i)=icol if(abs(a(icol,icol)).eq.0.0d0)then print*, 'Singular matrix ecountered in zgaussj [Case 2]' endif pivinv=1.0d0/a(icol,icol) a(icol,icol)=1.0d0 do 16 l=1,n a(icol,l)=a(icol,l)*pivinv 16 continue do 17 l=1,m b(icol,l)=b(icol,l)*pivinv 17 continue do 21 ll=1,n if(ll.ne.icol)then dum=a(ll,icol) a(ll,icol)=0. do 18 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 18 continue do 19 l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum 19 continue endif 21 continue 22 continue do 24 l=n,1,-1 if(indxr(l).ne.indxc(l))then do 23 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 23 continue endif 24 continue return end