Program CMTEST !Test package for CMTRIX complex*8 A0(40),A(40),B(40),C(40),D(40) complex*8 X0(10),W(10),X(10),Y(10),Z(10) complex*8 S integer*2 i,j,m,n,NF(10) 100 type 1000 1000 format('0','Test program for MATRIX package',/ 1 x,'Please enter m, n (2-6, m>n): ',$) accept 1005,m,n 1005 format(2I3) if(m.lt.2.or.m.gt.6)goto 100 if(n.lt.2.or.n.gt.6)goto 100 if(m*n.gt.40)goto 100 open(unit=1,name='CMTRIX.DAT',type='OLD') read(1,*) A0 read(1,*) X0 2000 format(5(2F8.4)) close(unit=1) type 1010 1010 format('0','MXWRT A: 1.0 to 12.0') call mxwrt(A0,m,n) type 1020 1020 format('0','MXMOV C: same as A') call mxmov(A0,m,n,C) call mxwrt(C,m,n) type 1030 1030 format('0','MXUNT I') call mxunt(B,m,n) call mxwrt(B,m,n) type 1040 1040 format('0','MXZER 0') call mxzer(B,m,n) call mxwrt(B,m,n) type 1050 1050 format('0','MXSCM B = 2*A') call mxmov(A0,m,n,B) s=cmplx(2.0,0.0) call mxscm(B,m,n,s) call mxwrt(B,m,n) type 1060 1060 format('0','MXADD D = A + C = 2*A') call mxadd(A0,C,m,n,D) call mxwrt(D,m,n) type 1070 1070 format('0','MXSUB C = D - B = 2*A-2*A = 0') call mxsub(D,B,m,n,C) call mxwrt(C,m,n) type 1080 1080 format('0','MXTRAP B = A''') call mxmov(A0,m,n,B) call mxtrap(B,m,n) call mxwrt(B,n,m) type 1090 1090 format('0','MXNORM') call mxnorm(A0,m,n,s) type 1100,s 1100 format(x,'Frobenius norm of A is (',F10.4,',',F10.4,')') type 1110 1110 format('0','MXMLT D = A * (A'')') call mxmov(A0,m,n,B) call mxtrap(B,m,n) call mxmlt(A0,B,m,n,m,D) call mxwrt(D,m,m) type 1120 1120 format('0','MXMBT C = A * A''') call mxmov(A0,m,n,B) call mxmbt(A0,B,m,n,m,C) call mxwrt(C,m,m) type 1130 1130 format('0','Check: B = D - C = 0') call mxsub(D,C,m,n,B) call mxwrt(B,m,m) type 1135 1135 format('0','MXMAT C = A'' * A') call mxmov(A0,m,n,B) call mxmat(B,A0,n,m,n,C) call mxwrt(C,n,n) type 1136 1136 format('0','Check: (A'')*A - A''*A should = 0') call mxtrap(B,m,n) call mxmlt(B,A0,n,m,n,D) call mxsub(D,C,n,n,B) call mxwrt(B,n,n) type 1140 1140 format('0','GINV A.X = C: find 1/A, & X') type 1150 1150 format(x,'Matrix A:') call mxmov(A0,m,n,A) call mxwrt(A,m,n) type 1160 1160 format('0','Vector C:') call mxmov(X0,n,1,Y) call mxwrt(Y,1,n) type * call ginv(A,m,n,D,Z,NF,Y,X,2,W) type 1170 1170 format('0','Inverse 1/A:') call mxwrt(A,n,m) type 1175 1175 format('0','Trimmed inverse 1/A:') call mxtrim(A,n,m,1.0e-6) call mxwrt(A,n,n) type 1180 1180 format('0','Solution vector X:') call mxwrt(X,n,1) type 1185 1185 format('0','Trimmed solution vector X:') call mxtrim(X,n,1,1.0e-6) call mxwrt(X,n,1) type 1190 1190 format('0','Check: B = 1/A * A, which should = I') call mxmlt(A,A0,n,m,n,B) call mxwrt(B,n,n) type 1200 1200 format('0','Trimmed B, which should = I exactly') call mxtrim(B,n,n,1.0e-5) call mxwrt(B,n,n) type 1210 1210 format('0','Check equation: A.X, which should = C') call mxmlt(A0,X,m,n,1,C) call mxwrt(C,m,1) type 1215 1215 format('0','Trimmed difference between A.X and C') call mxsub(X0,C,n,1,B) call mxtrim(B,n,1,1.0e-5) call mxwrt(B,n,1) goto 999 c This section not implemented: I don't understand the consequences! type 1220 1220 format('0','SVD: A = U.Q.V'' ') call mxmov(A0,m,n,A) type 1230 1230 format('0','Matrix A:') call mxwrt(A,m,n) type 1240 1240 format('0','Perform SVD itself') call SVD(A,m,n,B,C,D) type 1250 1250 format('0','Rank Orthonormal Matrix U') call mxwrt(B,m,m) type 1255 1255 format('0','Check: U.U'' =I') call mxmbt(B,B,m,m,m,A) call mxtrim(A,m,m,1.0e-5) call mxwrt(A,m,m) type 1270 1270 format('0','Orthogonal Matrix V') call mxwrt(C,m,m) type 1275 1275 format('0','Check: V.V'' =I') call mxmbt(C,C,n,n,n,A) call mxtrim(A,n,n,1.0e-5) call mxwrt(A,n,n) type 1280 1280 format('0','Matrix Q: Singular values of A''.A, ordered') call mxwrt(D,1,n) type 1290 1290 format('0','Verify A=U.Q.V'': calculate RHS') call mxzer(A,m,n) do 1292 i=1,min0(m,n) A((i-1)*n+i)=D(i) 1292 continue call mxmbt(A,C,m,m,m,D) call mxmlt(B,D,m,m,m,A) call mxtrim(A,m,m) call mxwrt(A,m,m) type 1300 1300 format('0','Trimmed difference A-U.Q.V'' ') call mxsub(A,A0,m,m,D) call mxtrim(D,m,m) call mxwrt(D,m,m) 999 type * stop end