C @(#)pca.for 17.1.1.1 (ES0-DMD) 01/25/02 17:16:39 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Carry out a PRINCIPAL COMPONENTS ANALYSIS C C (KARHUNEN-LOEVE EXPANSION). C C C C To call: CALL PCA(N,M,DATA,METHOD,IPRINT,A1,W1,W2,A2,IERR) where C C C C C C N, M : integer dimensions of ... C C DATA : input data. C C On output, DATA contains in first 7 columns the projections of C C the row-points on the first 7 principal components. C C METHOD: analysis option. C C = 1: on sums of squares & cross products matrix. C C = 2: on covariance matrix. C C = 3: on correlation matrix. C C IPRINT: print options. C C = 0: no printed output - arrays/vectors, only, contain items C C calculated. C C = 1: eigenvalues, only, output. C C = 2: printed output, in addition, of correlation (or other) C C matrix, eigenvalues and eigenvectors. C C = 3: full printing of items calculated. C C A1 : correlation, covariance or sums of squares & cross-products C C matrix, dimensions M * M. C C On output, A1 contains in the first 7 columns the projections C C of the column-points on the first 7 principal components. C C W1,W2 : real vectors of dimension M (see called routines for use). C C On output, W1 contains the cumulative percentage variances C C associated with the principal components. C C A2 : real array of dimensions M * M (see called routines for use). C C IERR : error indicator (normally zero). C C C C C C Inputs here are N, M, DATA, METHOD, IPRINT (and IERR). C C Output information is contained in DATA, A1, and W1. C C All printed outputs are carried out in easily recognizable subroutines C C called from the first subroutine following. C C C C Regarding precision and tolerances, we believe the eigen-extraction C C routine is reasonably precise; however, in subroutine PROJY we have C C taken a zero eigenvalue as being less than or equal to 0.00005 in C C value. This is adequate for most practical situations. C C C C C C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986. C C C C-------------------------------------------------------------------------C SUBROUTINE PCA(N,M,DATA,METHOD,IPRINT,A,W,FV1,Z,IERR) INTEGER N,M,METHOD,IPRINT,IERR,M2 REAL DATA(N,M), A(M,M), W(M), FV1(M), Z(M,M) C IF (METHOD.EQ.1) GOTO 100 IF (METHOD.EQ.2) GOTO 400 C If method.eq.3 or otherwise ... GOTO 700 C C Form sums of squares and cross-products matrix. C 100 CONTINUE CALL SCPCLP(N,M,DATA,A,IERR) IF (IERR.NE.0) GOTO 9000 C IF (IPRINT.GT.1) CALL OUTHMP(METHOD,M,A) C C Now do the PCA. C GOTO 1000 C C Form covariance matrix. C 400 CONTINUE CALL COVCLP(N,M,DATA,W,A,IERR) IF (IERR.NE.0) GOTO 9000 C IF (IPRINT.GT.1) CALL OUTHMP(METHOD,M,A) C C Now do the PCA. C GOTO 1000 C C Construct correlation matrix. C 700 CONTINUE CALL CORCLP(N,M,DATA,W,FV1,A,IERR) IF (IERR.NE.0) GOTO 9000 C IF (IPRINT.GT.1) CALL OUTHMP(METHOD,M,A) C C Now do the PCA. C GOTO 1000 C C Carry out eigenreduction. C 1000 M2 = M CALL TRED2P(M,M2,A,W,FV1,Z) CALL TQL2P(M,M2,W,FV1,Z,IERR) IF (IERR.NE.0) GOTO 9000 C C Output eigenvalues and eigenvectors. C IF (IPRINT.GT.0) CALL OUTVLP(N,M,W) IF (IPRINT.GT.1) CALL OUTVCP(N,M,Z) C C Determine projections and output them. C CALL PROJXP(N,M,DATA,Z,FV1) IF (IPRINT.EQ.3) CALL OUTPXP(N,M,DATA) CALL PROJYP(M,W,A,Z,FV1) IF (IPRINT.EQ.3) CALL OUTPYP(M,A) C 9000 RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine correlations of columns. C C First determine the means of columns, storing in WORK1. C C C C-------------------------------------------------------------------------C SUBROUTINE CORCLP(N,M,DATA,WORK1,WORK2,OUT,ISTAT) INTEGER N,M,ISTAT,I,J,J1,J2 REAL DATA(N,M), OUT(M,M), WORK1(M), WORK2(M) REAL EPS DATA EPS/1.E-10/ C DO 30 J = 1, M WORK1(J) = 0.0 DO 20 I = 1, N WORK1(J) = WORK1(J) + DATA(I,J) 20 CONTINUE WORK1(J) = WORK1(J)/FLOAT(N) 30 CONTINUE C C Next det. the std. devns. of cols., storing in WORK2. C DO 50 J = 1, M WORK2(J) = 0.0 DO 40 I = 1, N WORK2(J) = WORK2(J) + X (DATA(I,J)-WORK1(J))*(DATA(I,J)-WORK1(J)) 40 CONTINUE WORK2(J) = WORK2(J)/FLOAT(N) WORK2(J) = SQRT(WORK2(J)) IF (WORK2(J).LE.EPS) WORK2(J) = 1.0 50 CONTINUE C C Now centre and reduce the column points. C DO 70 I = 1, N DO 60 J = 1, M DATA(I,J) = (DATA(I,J)-WORK1(J))/(SQRT(FLOAT(N))*WORK2(J)) 60 CONTINUE 70 CONTINUE C C Finally calculate the cross product of the redefined data matrix. C DO 100 J1 = 1, M-1 OUT(J1,J1) = 1.0 DO 90 J2 = J1+1, M OUT(J1,J2) = 0.0 DO 80 I = 1, N OUT(J1,J2) = OUT(J1,J2) + DATA(I,J1)*DATA(I,J2) 80 CONTINUE OUT(J2,J1) = OUT(J1,J2) 90 CONTINUE 100 CONTINUE OUT(M,M) = 1.0 C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine covariances of columns. C C First determine the means of columns, storing in WORK. C C C C-------------------------------------------------------------------------C SUBROUTINE COVCLP(N,M,DATA,WORK,OUT,ISTAT) INTEGER N,M,ISTAT,I,J,J1,J2 REAL DATA(N,M), OUT(M,M), WORK(M) C DO 30 J = 1, M WORK(J) = 0.0 DO 20 I = 1, N WORK(J) = WORK(J) + DATA(I,J) 20 CONTINUE WORK(J) = WORK(J)/FLOAT(N) 30 CONTINUE C C Now centre the column points. C DO 50 I = 1, N DO 40 J = 1, M DATA(I,J) = DATA(I,J)-WORK(J) 40 CONTINUE 50 CONTINUE C C Finally calculate the cross product matrix of the redefined C data matrix. C DO 80 J1 = 1, M DO 70 J2 = J1, M OUT(J1,J2) = 0.0 DO 60 I = 1, N OUT(J1,J2) = OUT(J1,J2) + DATA(I,J1)*DATA(I,J2) 60 CONTINUE OUT(J2,J1) = OUT(J1,J2) 70 CONTINUE 80 CONTINUE C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine sums of squares and cross-products of columns. C C C C-------------------------------------------------------------------------C SUBROUTINE SCPCLP(N,M,DATA,OUT,ISTAT) INTEGER M,N,J1,J2,I,ISTAT REAL DATA(N,M), OUT(M,M) C DO 30 J1 = 1, M DO 20 J2 = J1, M OUT(J1,J2) = 0.0 DO 10 I = 1, N OUT(J1,J2) = OUT(J1,J2) + DATA(I,J1)*DATA(I,J2) 10 CONTINUE OUT(J2,J1) = OUT(J1,J2) 20 CONTINUE 30 CONTINUE C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Reduce a real, symmetric matrix to a symmetric, tridiagonal matrix. C C C C To call: CALL TRED2(NM,N,A,D,E,Z) where C C C C NM = row dimension of A and Z; C C N = order of matrix A (will always be <= NM); C C A = symmetric matrix of order N to be reduced to tridiagonal form; C C D = vector of dim. N containing, on output, diagonal elts. of trid. C C matrix. C C E = working vector of dim. at least N-1 to contain subdiagonal elts.; C C Z = matrix of dims. NM by N containing, on output, orthogonal C C transformation matrix producing the reduction. C C C C Normally a call to TQL2 will follow the call to TRED2 in order to C C produce all eigenvectors and eigenvalues of matrix A. C C C C Algorithm used: Martin et al., Num. Math. 11, 181-195, 1968. C C C C Reference: Smith et al., Matrix Eigensystem Routines - EISPACK C C Guide, Lecture Notes in Computer Science 6, Springer-Verlag, 1976, C C pp. 489-494. C C C C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986. C C C C-------------------------------------------------------------------------C SUBROUTINE TRED2P(NM,N,A,D,E,Z) INTEGER I,J,N,L,K,II,NM,JP1 REAL A(NM,N),D(N),E(N),Z(NM,N) REAL SCALE,F,G,H,HH C DO 100 I = 1, N DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE IF (N.EQ.1) GOTO 320 DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0 SCALE = 0.0 IF (L.LT.2) GOTO 130 DO 120 K = 1, L SCALE = SCALE + ABS(Z(I,K)) 120 CONTINUE IF (SCALE.NE.0.0) GOTO 140 130 E(I) = Z(I,L) GOTO 290 140 DO 150 K = 1, L Z(I,K) = Z(I,K)/SCALE H = H + Z(I,K)*Z(I,K) 150 CONTINUE C F = Z(I,L) G = -SIGN(SQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0 C DO 240 J = 1, L Z(J,I) = Z(I,J)/H G = 0.0 C Form element of A*U. DO 180 K = 1, J G = G + Z(J,K)*Z(I,K) 180 CONTINUE JP1 = J + 1 IF (L.LT.JP1) GOTO 220 DO 200 K = JP1, L G = G + Z(K,J)*Z(I,K) 200 CONTINUE C Form element of P where P = I - U U' / H . 220 E(J) = G/H F = F + E(J) * Z(I,J) 240 CONTINUE HH = F/(H + H) C Form reduced A. DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G DO 250 K = 1, J Z(J,K) = Z(J,K) - F*E(K) - G*Z(I,K) 250 CONTINUE 260 CONTINUE 290 D(I) = H 300 CONTINUE 320 D(1) = 0.0 E(1) = 0.0 C Accumulation of transformation matrices. DO 500 I = 1, N L = I - 1 IF (D(I).EQ.0.0) GOTO 380 DO 360 J = 1, L G = 0.0 DO 340 K = 1, L G = G + Z(I,K) * Z(K,J) 340 CONTINUE DO 350 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 350 CONTINUE 360 CONTINUE 380 D(I) = Z(I,I) Z(I,I) = 1.0 IF (L.LT.1) GOTO 500 DO 400 J = 1, L Z(I,J) = 0.0 Z(J,I) = 0.0 400 CONTINUE 500 CONTINUE C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine eigenvalues and eigenvectors of a symmetric, C C tridiagonal matrix. C C C C To call: CALL TQL2(NM,N,D,E,Z,IERR) where C C C C NM = row dimension of Z; C C N = order of matrix Z; C C D = vector of dim. N containing, on output, eigenvalues; C C E = working vector of dim. at least N-1; C C Z = matrix of dims. NM by N containing, on output, eigenvectors; C C IERR = error, normally 0, but 1 if no convergence. C C C C Normally the call to TQL2 will be preceded by a call to TRED2 in C C order to set up the tridiagonal matrix. C C C C Algorithm used: QL method of Bowdler et al., Num. Math. 11, C C 293-306, 1968. C C C C Reference: Smith et al., Matrix Eigensystem Routines - EISPACK C C Guide, Lecture Notes in Computer Science 6, Springer-Verlag, 1976, C C pp. 468-474. C C C C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986. C C C C-------------------------------------------------------------------------C SUBROUTINE TQL2P(NM,N,D,E,Z,IERR) INTEGER I,N,IERR,M,NM,L,J,L1,II,MML,K REAL D(N), E(N), Z(NM,N),EPS,F,B,H REAL P,R,G,C,S DATA EPS/1.E-12/ C IERR = 0 IF (N.EQ.1) GOTO 1001 DO 100 I = 2, N E(I-1) = E(I) 100 CONTINUE F = 0.0 B = 0.0 E(N) = 0.0 C DO 240 L = 1, N J = 0 H = EPS * (ABS(D(L)) + ABS(E(L))) IF (B.LT.H) B = H C Look for small sub-diagonal element. DO 110 M = L, N IF (ABS(E(M)).LE.B) GOTO 120 C E(N) is always 0, so there is no exit through the C bottom of the loop. 110 CONTINUE 120 IF (M.EQ.L) GOTO 220 130 IF (J.EQ.30) GOTO 1000 J = J + 1 C Form shift. L1 = L + 1 G = D(L) P = (D(L1)-G)/(2.0*E(L)) R = SQRT(P*P+1.0) D(L) = E(L)/(P+SIGN(R,P)) H = G-D(L) C DO 140 I = L1, N D(I) = D(I) - H 140 CONTINUE C F = F + H C QL transformation. P = D(M) C = 1.0 S = 0.0 MML = M - L C DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (ABS(P).LT.ABS(E(I))) GOTO 150 C = E(I)/P R = SQRT(C*C+1.0) E(I+1) = S * P * R S = C/R C = 1.0/R GOTO 160 150 C = P/E(I) R = SQRT(C*C+1.0) E(I+1) = S * E(I) * R S = 1.0/R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C Form vector. DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE 200 CONTINUE E(L) = S * P D(L) = C * P IF (ABS(E(L)).GT.B) GOTO 130 220 D(L) = D(L) + F 240 CONTINUE C C Order eigenvectors and eigenvalues. DO 300 II = 2, N I = II - 1 K = I P = D(I) DO 260 J = II, N IF (D(J).GE.P) GOTO 260 K = J P = D(J) 260 CONTINUE IF (K.EQ.I) GOTO 300 D(K) = D(I) D(I) = P DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE 300 CONTINUE C GOTO 1001 C Set error - no convergence to an eigenvalue after 30 iterations. 1000 IERR = 1 1001 RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output array. C C C C-------------------------------------------------------------------------C SUBROUTINE OUTMTP(N,M,ARRAY) INTEGER K1,K2,N,M REAL ARRAY(N,M) C DO 100 K1 = 1, N WRITE (6,1000) (ARRAY(K1,K2),K2=1,M) 100 CONTINUE C 1000 FORMAT(10(2X,F8.4)) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output half of (symmetric) array. C C C C-------------------------------------------------------------------------C SUBROUTINE OUTHMP(ITYPE,NDIM,ARRAY) INTEGER ITYPE,NDIM,K1,K2,STATUS REAL ARRAY(NDIM,NDIM) CHARACTER*80 MYLINE C IF (ITYPE.EQ.1) WRITE (MYLINE,1000) IF (ITYPE.EQ.2) WRITE (MYLINE,2000) IF (ITYPE.EQ.3) WRITE (MYLINE,3000) CALL STTPUT(MYLINE,STATUS) C DO 100 K1 = 1, NDIM WRITE (MYLINE,4000) (ARRAY(K1,K2),K2=1,K1) CALL STTPUT(MYLINE,STATUS) 100 CONTINUE C 1000 FORMAT(1X ,'SUMS OF SQUARES & CROSS-PRODUCTS MATRIX FOLLOWS.') 2000 FORMAT(1X ,'COVARIANCE MATRIX FOLLOWS.') 3000 FORMAT(1X ,'CORRELATION MATRIX FOLLOWS.') 4000 FORMAT(8(2X,F8.4)) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output eigenvalues in order of decreasing value. C C C C-------------------------------------------------------------------------C SUBROUTINE OUTVLP(N,NVALS,VALS) INTEGER K,NVALS,N,M,ISTAT REAL VALS(NVALS),TOT,CUM REAL VPC,VCPC CHARACTER LINE*80 C TOT = 0.0 DO 100 K = 1, NVALS TOT = TOT + VALS(K) 100 CONTINUE C C WRITE (6,1000) WRITE (LINE, 900) CALL STTPUT(LINE,ISTAT) WRITE(6,900) WRITE (LINE,1000) CALL STTPUT(LINE,ISTAT) WRITE(6,900) WRITE (LINE, 900) CALL STTPUT(LINE,ISTAT) CUM = 0.0 K = NVALS + 1 C M = NVALS C (We only want Min(nrows,ncols) eigenvalues output:) M = MIN0(N,NVALS) C C WRITE (6,1010) C WRITE (6,1020) WRITE (LINE,1010) CALL STTPUT(LINE,ISTAT) WRITE (LINE,1020) CALL STTPUT(LINE,ISTAT) 200 CONTINUE K = K - 1 CUM = CUM + VALS(K) VPC = VALS(K) * 100.0 / TOT VCPC = CUM * 100.0 / TOT C WRITE (6,1030) VALS(K),VPC,VCPC WRITE (LINE,1030) VALS(K),VPC,VCPC CALL STTPUT(LINE,ISTAT) C VALS(K) = VCPC IF (K.GT.NVALS-M+1) GOTO 200 C RETURN 900 FORMAT(' ') 1000 FORMAT(1X,' EIGENVALUES FOLLOW.') 1010 FORMAT X (' Eigenvalues As Percentages Cumul. Percentages') 1020 FORMAT X (' ----------- -------------- ------------------') 1030 FORMAT(F13.4,7X,F10.4,10X,F10.4) END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output FIRST SEVEN eigenvectors associated with eigenvalues C C in descending order. C C C C-------------------------------------------------------------------------C SUBROUTINE OUTVCP(N,NDIM,VECS) INTEGER NDIM,NUM,N,K1,K2 REAL VECS(NDIM,NDIM) C NUM = MIN0(N,NDIM,7) C WRITE(6,985) WRITE (6,1000) WRITE(6,985) WRITE (6,1010) WRITE (6,1020) DO 100 K1 = 1, NDIM WRITE (6,1030) K1,(VECS(K1,NDIM-K2+1),K2=1,NUM) 100 CONTINUE C RETURN 985 FORMAT(' ') 1000 FORMAT(1X ,'EIGENVECTORS FOLLOW.') 1010 FORMAT(' VBLE. EV-1 EV-2 EV-3 EV-4 EV-5 EV-6 X EV-7') 1020 FORMAT(' ------ ------ ------ ------ ------ ------ ------ X------') 1030 FORMAT(I5,2X,7F8.4) END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output projections of row-points on first 7 pricipal components. C C C C-------------------------------------------------------------------------C SUBROUTINE OUTPXP(N,M,PRJN) INTEGER N,M,K,NUM,J REAL PRJN(N,M) C NUM = MIN0(M,7) WRITE(6,985) WRITE (6,1000) WRITE(6,985) WRITE (6,1010) WRITE (6,1020) DO 100 K = 1, N WRITE (6,1030) K,(PRJN(K,J),J=1,NUM) 100 CONTINUE C 985 FORMAT(' ') 1000 FORMAT(1X ,'PROJECTIONS OF ROW-POINTS FOLLOW.') 1010 FORMAT(' OBJECT PROJ-1 PROJ-2 PROJ-3 PROJ-4 PROJ-5 PROJ-6 X PROJ-7') 1020 FORMAT(' ------ ------ ------ ------ ------ ------ ------ X ------') 1030 FORMAT(I5,2X,7F8.4) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output projections of columns on first 7 principal components. C C C C-------------------------------------------------------------------------C SUBROUTINE OUTPYP(M,PRJNS) INTEGER K,NUM,M,J REAL PRJNS(M,M) C NUM = MIN0(M,7) WRITE(6,985) WRITE (6,1000) WRITE(6,985) WRITE (6,1010) WRITE (6,1020) DO 100 K = 1, M WRITE (6,1030) K,(PRJNS(K,J),J=1,NUM) 100 CONTINUE C 985 FORMAT(' ') 1000 FORMAT(1X ,'PROJECTIONS OF COLUMN-POINTS FOLLOW.') 1010 FORMAT(' VBLE. PROJ-1 PROJ-2 PROJ-3 PROJ-4 PROJ-5 PROJ-6 X PROJ-7') 1020 FORMAT(' ------ ------ ------ ------ ------ ------ ------ X ------') 1030 FORMAT(I5,2X,7F8.4) RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Form projections of row-points on first 7 principal components. C C C C-------------------------------------------------------------------------C SUBROUTINE PROJXP(N,M,DATA,EVEC,VEC) INTEGER N,M,K,NUM,L,I,J REAL DATA(N,M), EVEC(M,M), VEC(M) C NUM = MIN0(M,7) DO 300 K = 1, N DO 50 L = 1, M VEC(L) = DATA(K,L) 50 CONTINUE DO 200 I = 1, NUM DATA(K,I) = 0.0 DO 100 J = 1, M DATA(K,I) = DATA(K,I) + VEC(J) * X EVEC(J,M-I+1) 100 CONTINUE 200 CONTINUE 300 CONTINUE C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine projections of column-points on 7 principal components. C C C C-------------------------------------------------------------------------C SUBROUTINE PROJYP(M,EVALS,A,Z,VEC) INTEGER M,NUM,L,J1,J2,J3 REAL EVALS(M), A(M,M), Z(M,M), VEC(M) C NUM = MIN0(M,7) DO 300 J1 = 1, M DO 50 L = 1, M VEC(L) = A(J1,L) 50 CONTINUE DO 200 J2 = 1, NUM A(J1,J2) = 0.0 DO 100 J3 = 1, M A(J1,J2) = A(J1,J2) + VEC(J3)*Z(J3,M-J2+1) 100 CONTINUE IF (EVALS(M-J2+1).GT.0.00005) A(J1,J2) = X A(J1,J2)/SQRT(EVALS(M-J2+1)) IF (EVALS(M-J2+1).LE.0.00005) A(J1,J2) = 0.0 200 CONTINUE 300 CONTINUE C RETURN END