C @(#)mda.for 17.1.1.1 (ES0-DMD) 01/25/02 17:16:38 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 MULTIPLE DISCRIMINANT ANALYSIS C C (DISCRIMINANT FACTOR ANALYSIS, C C CANONICAL DISCRININANT ANALYSIS). C C C C C C To call: CALL MDA(N,M,NG,DATA,GP,IPRINT,NOG,MEAN,MGP,TOTAL, C C BETWEEN,BETW2,CPROJ,W1,W2,IERR) where C C C C C C N, M : integer dimensions of ... C C DATA : input data (real). C C On output the first NG-1 columns of DATA contain projections C C of the N items on the discriminant factors. C C NG : (integer) number of groups. C C GP : Integer vector of length N giving group assignments. C C Must be specified correctly - no 0s or values > NG. C C IPRINT: integer; print options (= 3: full; otherwise none). C C NOG : integer vector of length NG (to contain group cardinalities).C C MEAN : real vector of length M (number of attributes or variables). C C MGP : real array of dimensions 2 by M. C C TOTAL : real array of dimensions M by M; on output contains inverse C C of total variance/covariance matrix. C C BETWEEN: real array of dimensions NG by NG. C C BETW2 : real array of dimensions NG by NG. C C CPROJ : real array of dimensions M by NG; on output contains the C C coefficients of the discriminant factors in terms of the C C original variables. C C W1, W2: real vectors of length M. C C IERR : initially 0; = 1 if there is no convergence in the TQL2 C C eigenroutine; = 2 if the total variance-covariance to be C C inverted is singular, - in this case, check that there are C C no columns with identical values, that N > M > NG, etc. C C C C C C Inputs here are N, M, NG, DATA, GP, IPRINT (and IERR). C C The principle output information is contained in DATA and CPROJ; and C C W1(NG-1), W1(NG-2) ... contain the eigenvalues in decreasing order. C C C C Notes: we require that N > M > NG (otherwise, an error is likely in C C the matrix inversion routine due to singularity). C C NG-1 eigenvalues eigenvectors are output. C C C C C C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986. C C C C-----------------------------------------------------------------------C SUBROUTINE MDA(N,M,NG,DATA,GP,IPRINT,NOG,MEAN,MGP,TOTAL, X BETWEEN,BETW2,CPROJ,W1,W2,IERR) REAL DATA(N,M),TOTAL(M,M),MEAN(M),MGP(NG,M),BETWEEN(NG,NG) REAL W1(M), W2(M), BETW2(NG,NG), CPROJ(M,NG) REAL D,D1,D2,D3,D4 INTEGER GP(N), NOG(NG), G,IMAT,IERR,NG2 INTEGER N,M,NG,IPRINT,I,J,J1,J2,K INTEGER K1,K2 C C Form global mean. C DO 200 J = 1, M MEAN(J) = 0.0 DO 100 I = 1, N MEAN(J) = MEAN(J) + DATA(I,J) 100 CONTINUE MEAN(J) = MEAN(J)/FLOAT(N) 200 CONTINUE C C Form (total) variance-covariance matrix. C DO 500 J1 = 1, M DO 400 J2 = 1, M TOTAL(J1,J2) = 0.0 DO 300 I = 1, N TOTAL(J1,J2) = TOTAL(J1,J2) + X (DATA(I,J1)-MEAN(J1))*(DATA(I,J2)-MEAN(J2)) 300 CONTINUE TOTAL(J1,J2) = TOTAL(J1,J2)/FLOAT(N) 400 CONTINUE 500 CONTINUE C IMAT = 1 C CALL OUTMTM(IMAT,M,TOTAL) C C Form group means. C DO 700 J = 1, M DO 600 K = 1, NG MGP(K,J) = 0.0 600 CONTINUE 700 CONTINUE C DO 900 I = 1, N G = GP(I) IF (G.EQ.0) GOTO 9000 NOG(G) = NOG(G) + 1 DO 800 J = 1, M MGP(G,J) = MGP(G,J) + DATA(I,J) 800 CONTINUE 900 CONTINUE C DO 1100 K = 1, NG DO 1000 J = 1, M MGP(K,J) = MGP(K,J)/NOG(K) 1000 CONTINUE 1100 CONTINUE C C Invert variance-covariance matrix. C CALL MINVM(M,TOTAL,D,W1,W2) IF (D.GT.0.000001) GOTO 1150 IERR = 2 GOTO 9000 1150 CONTINUE IMAT = 2 C CALL OUTMTM(IMAT,M,TOTAL) C C Form the symmetric variant of the between-groups C variance-covariance matrix for diagonalization. C DO 1200 K1 = 1, NG DO 1200 K2 = 1, NG BETWEEN(K1,K2) = 0.0 DO 1200 J1 = 1, M DO 1200 J2 = 1, M D1 = MGP(K1,J1) - MEAN(J1) D2 = FLOAT(NOG(K1))/FLOAT(N) D2 = SQRT(D2) D3 = MGP(K2,J2) - MEAN(J2) D4 = FLOAT(NOG(K2))/FLOAT(N) D4 = SQRT(D4) BETWEEN(K1,K2) = BETWEEN(K1,K2) + X (D1*D2)*TOTAL(J1,J2)*(D3*D4) 1200 CONTINUE C IMAT = 4 C CALL OUTMTM(IMAT,M,TOTAL) C C Carry out eigenreduction. C NG2 = NG CALL TRED2M(NG,NG2,BETWEEN,W1,W2,BETW2) CALL TQL2M(NG,NG2,W1,W2,BETW2,IERR) IF (IERR.NE.0) GOTO 9000 C C Output eigenvalues and eigenvectors. C IF (IPRINT.GT.1) CALL OUTEVM(N,M,NG,W1) IF (IPRINT.GT.1) CALL OUTVCM(N,M,NG,BETW2,NG-1) C C Convert eigenvectors in NG-space to those in M-space. C DO 1300 J = 1, M DO 1300 K = 1, NG CPROJ(J,K) = 0.0 DO 1300 J2 = 1, M DO 1300 K2 = 1, NG D1 = MGP(K2,J2) - MEAN(J2) D2 = FLOAT(NOG(K2))/FLOAT(N) D1 = D1*SQRT(D2) CPROJ(J,K)=CPROJ(J,K)+ X TOTAL(J,J2)*D1*BETW2(K2,NG-K+1) 1300 CONTINUE IF (IPRINT.GT.1) CALL OUTVCM(N,NG,M,CPROJ,NG-1) C C Determine projections and output them. C CALL PROJXM(N,M,NG,DATA,MEAN,CPROJ,W2,TOTAL) IF (IPRINT.EQ.3) CALL OUTPXM(N,M,NG,DATA) C C C 9000 CONTINUE RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE OUTMTM(IMAT,M,ARRAY) INTEGER M,K1,K2,IMAT REAL ARRAY(M,M) C WRITE (6,900) IMAT DO 100 K1 = 1, M WRITE (6,1000) (ARRAY(K1,K2),K2=1,M) 100 CONTINUE C 900 FORMAT(' IMAT =',I6) 1000 FORMAT(10(2X,F8.4)) RETURN END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Invert a symmetric matrix and calculate its determinant. C C C C C C To call: CALL MATINV(M,ARRAY,DET,W1,W2) where C C C C C C M : dimension of ... C C ARRAY : input matrix which is replaced by its inverse. C C NORDER : degree of matrix (order of determinant) C C DET : determinant of input matrix. C C W1, W2 : work vectors of dimension M. C C C C C C Reference: Philip B Bevington, "Data Reduction and Error Analysis C C for the Physical Sciences", McGraw-Hill, New York, 1969, C C pp. 300-303. C C C C C C F. Murtagh, ST-ECF, Garching-bei-Muenchen, January 1986. C C C C----------------------------------------------------------------------C SUBROUTINE MINVM(M,ARRAY,DET,IK,JK) INTEGER M,K,I,J,L REAL ARRAY(M,M), IK(M), JK(M) REAL DET,AMAX,SAVE C DET = 1.0 DO 100 K = 1, M C Find largest element ARRAY(I,J) in rest of matrix. AMAX = 0.0 21 DO 30 I = K, M DO 30 J = K, M IF (ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30 24 AMAX = ARRAY(I,J) IK(K) = I JK(K) = J 30 CONTINUE C Interchange rows and columns to put AMAX in ARRAY(K,K). IF (AMAX) 41,32,41 32 DET = 0.0 GOTO 140 41 I = IK(K) IF (I-K) 21,51,43 43 DO 50 J = 1, M SAVE = ARRAY(K,J) ARRAY(K,J) = ARRAY(I,J) 50 ARRAY(I,J) = -SAVE 51 J = JK(K) IF (J-K) 21,61,53 53 DO 60 I = 1, M SAVE = ARRAY(I,K) ARRAY(I,K) = ARRAY(I,J) 60 ARRAY(I,J) = -SAVE C Accumulate elements of inverse matrix. 61 DO 70 I = 1, M IF (I-K) 63,70,63 63 ARRAY(I,K) = -ARRAY(I,K)/AMAX 70 CONTINUE DO 80 I = 1, M DO 80 J = 1, M IF (I-K) 74,80,74 74 IF (J-K) 75,80,75 75 ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K)*ARRAY(K,J) 80 CONTINUE DO 90 J = 1, M IF (J-K) 83,90,83 83 ARRAY(K,J) = ARRAY(K,J)/AMAX 90 CONTINUE ARRAY(K,K) = 1.0/AMAX 100 DET = DET * AMAX C Restore ordering of matrix. DO 130 L = 1, M K = M - L + 1 J = IK(K) IF (J-K) 111,111,105 105 DO 110 I = 1, M SAVE = ARRAY(I,K) ARRAY(I,K) = -ARRAY(I,J) 110 ARRAY(I,J) = SAVE 111 I = JK(K) IF (I-K) 130,130,113 113 DO 120 J = 1, M SAVE = ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J) = SAVE 130 CONTINUE 140 RETURN END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Reduce a real, symmetric matrix to a symmetric, tridiagonal C C 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. 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 producting 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 Garching-bei-Muenchen, January 1986. C C C C----------------------------------------------------------------------------C SUBROUTINE TRED2M(NM,N,A,D,E,Z) C INTEGER N,NM,I,J,II,L,K,JP1,HH REAL A(NM,N),D(N),E(N),Z(NM,N) REAL SCALE,F,G,H 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 Garching-bei-Muenchen, January 1986. C C C C----------------------------------------------------------------------------C SUBROUTINE TQL2M(NM,N,D,E,Z,IERR) C INTEGER N,NM,IERR,I,L,J,M,L1,MML,II,K REAL D(N), E(N), Z(NM,N) REAL EPS,F,B,H,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 bottom C 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 eigenvalues in order of decreasing value. C C C C-----------------------------------------------------------------------C SUBROUTINE OUTEVM(N,M,NG,VALS) INTEGER N,M,NG,K REAL VALS(NG),TOT,CUM,VPC,VCPC C TOT = 0.0 DO 100 K = 2, NG TOT = TOT + VALS(K) 100 CONTINUE C WRITE (6,1000) CUM = 0.0 K = NG + 1 WRITE (6,1010) WRITE (6,1020) 200 CONTINUE K = K - 1 CUM = CUM + VALS(K) VPC = VALS(K) * 100.0 / TOT VCPC = CUM * 100.0 / TOT WRITE (6,1030) VALS(K),VPC,VCPC IF (K.GT.2) GOTO 200 C RETURN 1000 FORMAT(1H0,'EIGENVALUES FOLLOW.',/) 1010 FORMAT(' Eigenvalues As Percentages Cumul. Percentages') 1020 FORMAT(' ----------- -------------- ------------------') 1030 FORMAT(F10.4,9X,F10.4,10X,F10.4) END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Output FIRST SEVEN eigenvectors associated with eigenvalues in C C decreasing order. C C C C-----------------------------------------------------------------------C SUBROUTINE OUTVCM(N1,N2,N3,VECS,N4) INTEGER N1,N2,N3,N4,NUM,K1,K2 REAL VECS(N3,N3) C NUM = MIN0(N4,7) WRITE (6,1000) WRITE (6,1010) WRITE (6,1020) DO 100 K1 = 1, N3 WRITE (6,1030) K1,(VECS(K1,K2),K2=1,NUM) 100 CONTINUE C RETURN 1000 FORMAT(1H0,'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 on discriminant factors. C C C C-----------------------------------------------------------------------C SUBROUTINE OUTPXM(N,M,NG,PRJN) INTEGER N,M,NUM,NG,K,J REAL PRJN(N,M) C NUM = MIN0(N,M,NG,7) WRITE (6,1000) WRITE (6,1010) WRITE (6,1020) DO 100 K = 1, N WRITE (6,1030) K,(PRJN(K,J),J=1,NUM-1) 100 CONTINUE C 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 column points on up to first 7 discriminant C C axes. C C C C-----------------------------------------------------------------------C SUBROUTINE OUTPYM(N,M,NG,PRJNS) INTEGER N,M,NG,K,NUM,J REAL PRJNS(NG,NG) C NUM = MIN0(N,M,NG,7) WRITE (6,1000) WRITE (6,1010) WRITE (6,1020) DO 100 K = 1, M WRITE (6,1030) K,(PRJNS(K,J),J=1,NUM) 100 CONTINUE C 1000 FORMAT(1H0,'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 (up to) first 7 factors. C C C C-----------------------------------------------------------------------C SUBROUTINE PROJXM(N,M,NG,DATA,MEAN,EVEC,VEC,TOTINV) INTEGER N,M,K,L,I,J1,NG,NUM REAL DATA(N,M), EVEC(M,M), VEC(M), TOTINV(M,M), MEAN(M) C NUM = MIN0(N,M,NG,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 J1 = 1, M CC DO 75 J2 = 1, M DATA(K,I) = DATA(K,I) + (VEC(J1) - MEAN(J1))* X EVEC(J1,I) CC 75 CONTINUE 100 CONTINUE 200 CONTINUE 300 CONTINUE C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C C C C Determine projections of column points on (up to) 7 factors. C C C C-----------------------------------------------------------------------C SUBROUTINE PROJYM(N,M,NG,EVALS,A,Z,VEC) INTEGER M,N,NG,NUM,L,J1,J2,J3 REAL EVALS(M), A(M,M), Z(M,M), VEC(M) C NUM = MIN0(N,M,NG,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.0) A(J1,J2) = X A(J1,J2)/SQRT(EVALS(M-J2+1)) IF (EVALS(M-J2+1).EQ.0.0) A(J1,J2) = 0.0 200 CONTINUE 300 CONTINUE C RETURN END