c CMTRIX: A package for matrix manipulation for complex variables, c and especially for the least squares solution of over-determined but c noisy linear equations. c c It demonstrates one of the beauties of Fortran: arrays of unknown c dimensions in mathematical libraries. c c Real version (MATRIX) founded on HWH's Matrix package (11-471) c c Complex version derived from the above by R N Caffin Aug-86 c c c Subroutine MXZER(A,m,n) c sets A(m,n) up as a zero matrix. complex*8 a(1) integer*2 i,m,n do 10 i = 1,m*n a(i) = cmplx(0.0,0.0) 10 continue return end Subroutine MXUNT(A,m,n) c sets matrix A(m,n) up with ones in the leading diagonal and c with zeroes elsewhere. complex*8 a(1) integer*2 i,m,n do 10 i = 1,m*n a(i) = cmplx(0.0,0.0) !Zero all 10 continue do 20 i = 1,min0(m,n) a(i-n+n*i ) = cmplx(1.0,0.0) !Then set diagonal 20 continue return end Subroutine MXSCM(A,m,n,s) c multiplies the matrix A(m,n) by the complex scalar factor s . complex*8 s,a(1) integer*2 i,m,n do 10 i = 1,m*n a(i) = s*a(i) 10 continue return end Subroutine MXADD(A,B,m,n,R) c adds matrix B(m,n) to matrix A(m,n) and returns the result as c R(m,n). If desired R in the call may be replaced by A or B,then c the result overwrites A or B as the case may be. complex*8 a(1),b(1),r(1) integer*2 i,m,n do 10 i = 1,m*n r(i) = a(i) + b(i) 10 continue return end Subroutine MXSUB(A,B,m,n,R) c subtracts the matrix B(m,n) from the matrix A(m,n) and returns c the result as R(m,n). If desired R in the call may be replaced c by A or B ,then the result overwrites A or B as the case may be. complex*8 a(1),b(1),r(1) integer*2 i,m,n do 10 i = 1,m*n r(i) = a(i) - b(i) 10 continue return end Subroutine MXMOV(A,m,n,R) c transfers matrix A(m,n) into a single array R( ) row by row. complex*8 a(1),r(1) integer*2 i,m,n do 10 i=1,m*n r(i) = a(i) 10 continue return end Subroutine MXNORM(A,m,n,s) c Calculates the (complex?) Frobenius norm of a(m,n) and returns it c as complex s. complex*8 a(1) complex*8 dot,s integer*2 i,m,n s = cmplx(0.0,0.0) do 10 i = 1,n s = s + dot(A,m,n,i,i) 10 continue s = sqrt(s) return end Subroutine MXMLT(A,B,m,l,n,R) c Multiply A(m,l) by B(l,n) and store result of order (m,n) in R( ) complex*8 a(1),b(1),r(2) integer*2 l,m,n r(1) = (0.0,0.0) r(2) = (0.0,0.0) call mlt(a,b,m,l,n,r) return end Subroutine MXMBT(A,B,m,l,n,R) c Multiply A(m,l) by transpose of B(n,l) to give product of order c (m,n) to be stored in R( ) . complex*8 a(1),b(1),r(2) integer*2 l,m,n r(1) = (0.0,0.0) r(2) = (1.0,0.0) call mlt(a,b,m,l,n,r) return end Subroutine MXMAT(A,B,m,l,n,R) c Multiplies the product of the transpose of A(l,m) with B(l,n) and c stores the product of order (m,n) in R( ). complex*8 a(1),b(1),r(2) integer*2 l,m,n r(1) = (1.0,0.0) r(2) = (0.0,0.0) call mlt(a,b,m,l,n,r) return end Subroutine MLT(A,B,m,l,n,R) c Internal routine called by MXMLT, MXMBT, MXMAT complex*8 a(1),b(1),r(2) integer*2 i,j,k,l,m,n integer*2 k1,k2,kml,kln,kr,kki,kjk,ka,kb k1 = ifix(real(r(1))) k2 = ifix(real(r(2))) kml=k1*(m-l) kln=k2*(l-n) do 120 i=1,m do 110 j=1,n kr=(i-1)*n+j r(kr)=cmplx(0.0,0.0) do 100 k=1,l kki=k1*(k-i) kjk=k2*(j-k) ka=(kki+i-1)*(kml+l)-kki+k kb=(kjk+k-1)*(kln+n)-kjk+j r(kr)=r(kr)+a(ka)*b(kb) 100 continue 110 continue 120 continue return end Subroutine MXTRIM(r,m,n,eps) c Sets to zero those elements of R smaller in absolute value than eps complex*8 r(1) real*4 a,b,eps integer*2 i,m,n do 100 i = 1,m*n a=real(r(i)) b=aimag(r(i)) if(abs(a).lt.eps) a = 0.0 if(abs(b).lt.eps) b = 0.0 r(i)=cmplx(a,b) 100 continue return end Subroutine MXWRT(r,m,n) c Types out elements of a one row to a line up to 6 elements per c line,with spill over to next line for more than 6 elements per c line. complex*8 r(1) integer*2 i,j,m,n integer*2 k type 1000 1000 format(x) do 100 i=1,m k=(i-1)*n type 1010, i,((real(r(j)),aimag(r(j))),j=k+1,k+min0(3,n)) 1010 format(x,I2,3(x,'(',G11.5,',',G11.5,')',:)) if(n.gt.3)type 1020, ((real(r(j)),aimag(r(j))),j=k+4,k+n) 1020 format(x,2x,3(x,'(',G11.5,',',G11.5,')',:)) 100 continue return end Subroutine MXPRT(A,m,n) c Prints out elements of a one row to a line up to 5 (complex) elements c per line ,with spill over to next line for more than 5 elements per c line. complex*8 A(1) integer*2 i,j,m,n integer*2 k print 1000 1000 format(x) do 100 i=1,m k=(i-1)*n print 1010, i,((real(a(j)),aimag(a(j))),j=k+1,k+min0(5,n)) 1010 format(x,I2,5(x,'(',G11.5,',',G11.5,')',)) if(n.gt.5)print 1020, ((real(a(j)),aimag(a(j))),j=k+6,k+n) 1020 format(x,2x,5(x,'(',G11.5,',',G11.5,')',)) 100 continue return end Subroutine MXTRAP(A,m,n) c Returns the transpose of matrix A(m,n) stored row by row in a c single row array A as A'(n,m) stored similarly in array A. complex*8 sav,A(1) mn = m*n-1 mnn = mn-n mn1 = mn-1 k = 0 i1 = 1 do 140 i = 2,mn1 if(k.gt.mnn) go to 100 k = k+n j = k go to 110 100 k = k-mnn j = k 110 if(j.ge.i1) go to 120 jn = j*n j = jn-mn*(jn/mn) go to 110 120 if(j.eq.i1) go to 130 j = j+1 sav = a(i) a(i) = a(j) a(j) = sav 130 i1 = i 140 continue return end complex function DOT(R,m,n,j,k) c Computes the inner product of columns j and k in double precision c and returns the value as function dot in single precision . complex*8 R(1) real*8 a,b,c,d,sr,si sr=0.0 si=0.0 do 100 i = 1,m ij = (i-1)*n+j ik = ij+k-j a=real(R(ij)) !complex*8 to real*4 to real*8 b=aimag(R(ij)) c=real(R(ik)) d=aimag(R(ik)) sr=sr+a*c-b*d !Sum real part in DP si=si+a*d+b*c !Sum imaginary part in DP 100 continue dot = cmplx(sngl(sr),sngl(si)) return end complex function DOT1(U,m,n,V,j,k,m0) complex*8 U(1),V(1) real*8 a,b,c,d,sr,si sr=0.0 si=0.0 if(m0.gt.m) goto 110 do 100 i = m0,m ij = (i-1)*n+j ik = ij+k-j a=real(U(ij)) !complex*8 to real*4 to real*8 b=aimag(U(ij)) c=real(V(ik)) d=aimag(V(ik)) sr=sr+a*c-b*d !Sum real part in DP si=si+a*d+b*c !Sum imaginary part in DP 100 continue 110 dot1 = cmplx(sngl(sr),sngl(si)) return end complex function DOT2(U,m,n,V,j,k,n0) complex*8 U(1),V(1) real*8 a,b,c,d,sr,si sr=0.0 si=0.0 if(n0.gt.n) goto 110 do 100 i = n0,n ji = -n+n*j+i ki = -n+n*k+i a=real(U(ji)) !complex*8 to real*4 to real*8 b=aimag(U(ji)) c=real(V(ki)) d=aimag(V(ki)) sr=sr+a*c-b*d !Sum real part in DP si=si+a*d+b*c !Sum imaginary part in DP 100 continue 110 dot2 = cmplx(sngl(sr),sngl(si)) return end complex function DOT3(U,m,n,V,i,j,n0) complex*8 U(1),V(1) real*8 a,b,c,d,sr,si sr=0.0 si=0.0 if(n0.gt.n) return do 100 k = n0,n ik = -n+n*i+k kj = -n+j+n*k a=real(U(ik)) !complex*8 to real*4 to real*8 b=aimag(U(ik)) c=real(V(kj)) d=aimag(V(kj)) sr=sr+a*c-b*d !Sum real part in DP si=si+a*d+b*c !Sum imaginary part in DP 100 continue 110 dot3 = cmplx(sngl(sr),sngl(si)) return end Subroutine GINV(A,m,n,U,AT,NF,C,X,ix,CQ) c c This Subroutine calculates the (Moore-Penrose) Generalized Inverse c or Pseudo Inverse A+(n,m) of the matrix A(m,n), m>=n . c It solves A.x=C either directly (m=n) or in a least squares sense. c See (eg) Stewart G W, Introduction to Matrix Computations. c c The result is stored in A(m*n) space. c m is the number of rows in A(m,n) . c n is the number of columns in A(m,n) . c U(n,n) is the book-keeping matrix. c ix is to be set by the user depending upon the intended use of the c algorithm,see below. c If(ix = 0,1,2) computes A+(n,m) or X(n) or both respectively. c A,U,AT,NF,C,X,CQ are all singly dimensioned arrays which are c appropriately dimensioned in the calling program. c A the coefficient array is of dimension m*n and is entered row by c row in the main program. c U the book-keeping array is of dimension n*n . c U,AT,NF, and CQ are included as formal parameters only to permit c their dimensions to be set in the main program. AT and NF are of c dimension n-1 and CQ is of dimension n . c The vector of constant coefficients C is of dimension m and is c supplied by the main program. If ix is less than 1, C need not be c given by the main program and may be of dimension 1 . c The vector of unknowns X is of dimension n, but may be of c dimension 1,only, if ix is less than 1 . c complex*8 A(1),U(1),AT(1),C(1),X(1),CQ(1) complex*8 dot,dota,dotb,fac,sav complex*8 one,zero real*4 hf,th,eps,tol,far integer*2 NF(1) one = cmplx(1.0,0.0) zero = cmplx(0.0,0.0) hf = 0.5 th = 200.0 kr = 0 c c dependent column tolerance for nb bit floating fraction . c nb = 24 eps = hf**(nb+1) tol = th*eps tol = tol*tol do 230 j = 1,n jm1 = j-1 dota = dot(a,m,n,j,j) nf(j) = 1 do 110 i = 1,n ij = (i-1)*n+j u(ij) = zero !Enter latest linearized a(ij) values. 110 continue jj = (j-1)*n+j u(jj) = one if(jm1.eq.0) go to 160 c c Gram-Schmidt orthogonalization with correction if required . c fac = (eps/float(m))*sqrt(dota) far = cabs(fac) ll = 0 120 ll = ll+1 lt = 1 do 130 k = 1,jm1 at(k) = dot(a,m,n,j,k) 130 continue do 150 k = 1,jm1 if(cabs(at(k)).gt.far) lt = ll-3 do 140 i = 1,m ij = (i-1)*n+j ik = ij+k-j if(nf(k).ne.0) a(ij) = a(ij) - at(k)*a(ik) if(i.le.k) u(ij) = u(ij) -at(k)*u(ik) 140 continue 150 continue if(lt.le.0) go to 120 160 dotb = dot(a,m,n,j,j) if(cabs(dota).lt.1.0e-37) go to 190 if(cabs(dotb/dota).gt.tol) go to 190 c c modified treatment where column is not effectively independent . c do 170 k = 1,jm1 at(k) = dot(u,k,n,j,k) 170 continue do 180 i = 1,m ij = (i-1)*n+j do 180 k = 1,jm1 if(nf(k).eq.0) go to 180 if(cabs(at(k)).le.eps) go to 180 ik = ij+k-j a(ij) = a(ij) - at(k)*a(ik) 180 continue nf(j) = 0 fac = dot(u,j,n,j,j) fac = one/sqrt(fac) go to 200 c c kr gives the total of independent columns, i.e. column rank at end c 190 fac = one/sqrt(dotb) kr = kr+1 c c column normalized. c 200 do 210 i = 1,m ij = (i-1)*n+j a(ij) = fac*a(ij) if(i.le.j) u(ij) = fac*u(ij) 210 continue if(ix.lt.1) go to 230 fac = zero do 220 i = 1,m ij = (i-1)*n+j fac = fac + c(i)*a(ij) !Enter latest linearized c(i) values. 220 continue cq(j) = fac 230 continue if(kr.lt.n) type 1000, kr,n 1000 format(x,'Matrix A(m,n) is of rank',I4,' with',i3,' columns') if(m.eq.n) go to 240 go to 260 240 sav = one do 250 i = 1,n ii = (i-1)*n+i sav = sav*u(ii) 250 continue sav = one/sav type 1010, real(sav),aimag(sav) 1010 format(x,'Value of determinant is ',G16.8,x,G16.8) 260 if(ix.lt.0) return if(ix.eq.0) go to 290 do 280 j = 1,n fac = zero do 270 k = j,n jk = (j-1)*n+k fac = fac + cq(k)*u(jk) 270 continue x(j) = fac 280 continue if(ix.eq.1) return 290 do 310 j = 1,n do 310 i = 1,m ij = (i-1)*n+j fac = zero do 300 k = j,n ik = ij+k-j jk = (j-1)*n+k fac = fac + a(ik)*u(jk) 300 continue a(ij) = fac 310 continue call mxtrap(a,m,n) return end Subroutine SVD(A,m,n,U,V,Q) c c Currently NOT supported c c This computes the singular values and complete orthogonal c decomposition of the real matrix A(m,n) to A = U*Q*V', where c U'*U = V'*V = I, m >= n. c The number of non-zero elements in Q gives the rank of A. c For a positive definite matrix the singular values are the same as c the eigenvalues. c See (eg) Stewart G W, Introduction to Matrix Computations. c c Matrix Q contains the singular values of (A'*A) so ordered that c q(i) >= q(i+1), i = 1,n-1. Q() is the diagonal matrix of c singular values. Eps specifies the machine precision, that is c eps is the smallest value for which (1.0+eps) > 1.0 in the c computer being used. Tol >= amin/eps, where amin is the c smallest positive number for which amin > 0.0 . c complex*8 a(1),u(1),v(1),q(1),e(81) complex*8 c,f,g,h,s,x,y,z complex*8 z0,w complex*8 dot,dot1,dot2,dot3 type *,'Complex SVD not supported' call mxzer(U,m,m) call mxzer(V,m,m) call mxzer(Q,m,m) return c********************************************************************** n1 = n+1 z0 = cmplx(0.0,0.0) hf = 0.5 w = cmplx(1.0,0.0) amin = 1.0e-38 nb = 24 eps = hf**(nb) tol = amin/eps mn = m*n do 5 i = 1,mn u(i) = a(i) 5 continue c c Householder's reduction to bidiagonal form. c x = z0 g = z0 do 130 i = 1,n e(i) = g l = i+1 ii = (-n+i)+n*i s = dot1(u,m,n,u,i,i,i) if(abs(s).ge.tol) go to 10 g = z0 go to 50 10 f = u(ii) g = sqrt(s) !?? c*** if(f.lt.z0) go to 20 !???? g = -g 20 h = f*g-s u(ii) = f-g if(l.gt.n) go to 50 do 40 j = l,n s = dot1(u,m,n,u,i,j,i) f = s/h do 30 k = i,m ki = (i-n)+n*k kj = (-i)+ki+j u(kj) = u(kj)+f*u(ki) 30 continue 40 continue 50 continue q(i) = g s = dot2(u,m,n,u,i,i,l) if(abs(s).ge.tol) go to 60 g = z0 go to 120 60 f = u(ii+1) g = sqrt(s) !?? c*** if(f.lt.z0) go to 70 !?????? g = -g 70 h = f*g-s u(ii+1) = f-g if(l.gt.n) go to 90 do 80 j = l,n ij = ii-i+j e(j) = u(ij)/h 80 continue 90 if(l.gt.m) go to 120 do 110 j = l,m s = dot2(u,m,n,u,j,i,l) do 100 k = l,n jk = (-n+n*j)+k u(jk) = u(jk)+s*e(k) 100 continue 110 continue 120 continue y = abs(q(i)) + abs(e(i)) !??????? c**** if(y.gt.x) x = y !??????? 130 continue c c accumulation of right hand transformations. c do 200 in = 1,n i = n1-in l = i+1 ii = (-n+i)+n*i c*** if(g.eq.z0) go to 170 !??? h = g*u(ii+1) if(l.gt.n) go to 190 do 140 j = l,n ji = (i-n)+n*j ij = -n+n*i+j v(ji) = u(ij)/h 140 continue do 160 j = l,n s = dot3(u,m,n,v,i,j,l) do 150 k = l,n ki = (-n+i)+n*k kj = (j-i)+ki v(kj) = v(kj)+s*v(ki) 150 continue 160 continue 170 continue do 180 j = l,n ij = (-n+n*i)+j ji = (-n+i)+n*j v(ji) = z0 v(ij) = z0 180 continue 190 continue v(ii) = w g = e(i) 200 continue c c accumulation of left hand transformations. c do 300 in = 1,n i = n1-in l = i+1 ii = -n+n*i+i g = q(i) if(l.gt.n) go to 220 do 210 j = l,n ij = ii-i+j u(ij) = z0 210 continue 220 continue c**** if(g.eq.z0) go to 270 !????? h = g*u(ii) if(l.gt.n) go to 250 do 240 j = l,n s = dot1(u,m,n,u,i,j,l) f = s/h do 230 k = i,m ki = -n+i+n*k kj = j-i+ki u(kj) = u(kj)+f*u(ki) 230 continue 240 continue 250 continue do 260 j = i,m ji = -n+i+n*j u(ji) = u(ji)/g 260 continue go to 290 270 do 280 j = i,m ji = -n+i+n*j u(ji) = z0 280 continue 290 u(ii) = u(ii)+w 300 continue c c diagonalization of the bidiagonal form. c eps = abs(x)*eps !???? orig x*eps do 440 kn = 1,n k = n1-kn 310 do 320 ln = 1,k !Test for f splitting. l = 1+k-ln if(abs(e(l)).le.eps) go to 360 !?? if(l.lt.2) go to 320 if(abs(q(l-1)).le.eps) go to 330 !?? 320 continue 330 c = z0 !Cancellation of e(l) if l > 1 s = w l1 = l-1 if(l.lt.1.or.l.gt.k) go to 360 do 350 i = l,k f = s*e(i) e(i) = c*e(i) if(abs(f).le.eps) go to 360 !?? g = q(i) h = sqrt(f*f+g*g) !??? q(i) = h c = g/h s = -f/h do 340 j = 1,m ji = -n+i+n*j jl1 = l1-i+ji y = u(jl1) z = u(ji) u(jl1) = y*c+z*s u(ji) = -y*s+z*c 340 continue 350 continue !Cancellation of e(l) if l > 1 360 z = q(k) !Test for convergence. if(l.eq.k) go to 410 x = q(l) !Shift from bottom (2*2) minor. y = q(k-1) g = e(k-1) h = e(k) f = ((y-z)*(y+z)+(g-h)*(g+h))/(h*(y+y)) g = sqrt(f*f+w) !???? gs = g c**** if(f.lt.z0) gs = -g !????? f = ((x-z)*(x+z)+h*(y/(f+gs)-h))/x s = w !next qr transformation. c = w lp1 = l+1 if(lp1.gt.k) go to 400 do 390 i = lp1,k g = e(i) y = q(i) h = s*g g = c*g z = sqrt(f*f+h*h) !??? e(i-1) = z c = f/z s = h/z f = x*c+g*s g = -x*s+g*c h = y*s y = y*c do 370 j = 1,n ji = -n+i+n*j x = v(ji-1) z = v(ji) v(ji-1) = x*c+z*s v(ji) = -x*s+z*c 370 continue z = sqrt(f*f+h*h) !??? q(i-1) = z c = f/z s = h/z f = c*g+s*y x = -s*g+c*y do 380 j = 1,m ji = -n+i+n*j y = u(ji-1) z = u(ji) u(ji-1) = y*c+z*s u(ji) = -y*s+z*c 380 continue e(l) = z0 390 continue 400 continue e(k) = f q(k) = x go to 310 c c convergence! c 410 continue c**** if(z.ge.z0) go to 430 c c make q(k) non-negative. c q(k) = -z do 420 j = 1,n jk = -n+k+n*j v(jk) = -v(jk) 420 continue 430 continue 440 continue n1 = n-1 c c sort on the singular values of (a'.a) and related eigenvectors. c do 490 k = 1,n1 g = -w j = k do 450 i = k,n c**** if(q(i).le.g) go to 450 !???? g = q(i) j = i 450 continue if(j.eq.k) go to 480 q(j) = q(k) q(k) = g do 460 i = 1,n ij = -n+n*i+j ik = ij+k-j s = v(ij) v(ij) = v(ik) v(ik) = s 460 continue do 470 i = 1,m ij = -n+n*i+j ik = ij+k-j s = u(ij) u(ij) = u(ik) u(ik) = s 470 continue 480 continue 490 continue c c Zero the effectively zero singular values. c call mxtrim(q,1,n,eps) return end