* SUBROUTINE MXDCMM ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A * BY A VECTOR X. * * PARAMETERS : * II N NUMBER OF ROWS OF THE MATRIX A. * II M NUMBER OF COLUMNS OF THE MATRIX A. * RI A(N*M) RECTANGULAR MATRIX STORED COLUMNWISE IN THE * ONE-DIMENSIONAL ARRAY. * RI X(M) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR EQUAL TO A*X. * * SUBPROGRAMS USED : * S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. * S MXVSET INITIATION OF A VECTOR. * SUBROUTINE MXDCMM(N,M,A,X,Y) INTEGER N,M DOUBLE PRECISION A(*),X(*),Y(*) INTEGER J,K CALL MXVSET(N,0.0D 0,Y) K=0 DO 1 J=1,M CALL MXVDIR(N,X(J),A(K+1),Y,Y) K=K+N 1 CONTINUE RETURN END * SUBROUTINE MXDPGB ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC * POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L) * OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RU X(N) ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR * EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR * EQUATIONS. * II JOB OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN * X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X. * * METHOD : * BACK SUBSTITUTION * SUBROUTINE MXDPGB(N,A,X,JOB) INTEGER JOB,N DOUBLE PRECISION A(*),X(*) INTEGER I,II,IJ,J IF (JOB.GE.0) THEN * * PHASE 1 : X:=L**(-1)*X * IJ = 0 DO 20 I = 1,N DO 10 J = 1,I - 1 IJ = IJ + 1 X(I) = X(I) - A(IJ)*X(J) 10 CONTINUE IJ = IJ + 1 20 CONTINUE END IF IF (JOB.EQ.0) THEN * * PHASE 2 : X:=D**(-1)*X * II = 0 DO 30 I = 1,N II = II + I X(I) = X(I)/A(II) 30 CONTINUE END IF IF (JOB.LE.0) THEN * * PHASE 3 : X:=TRANS(L)**(-1)*X * II = N* (N-1)/2 DO 50 I = N - 1,1,-1 IJ = II DO 40 J = I + 1,N IJ = IJ + J - 1 X(I) = X(I) - A(IJ)*X(J) 40 CONTINUE II = II - I 50 CONTINUE END IF RETURN END * SUBROUTINE MXDPGD ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF A DIRECTION OF NEGATIVE CURVATURE WITH RESPECT TO A * DENSE SYMMETRIC MATRIX A USING THE FACTORIZATION A+E=L*D*TRANS(L) * OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RO X(N) COMPUTED DIRECTION OF NEGATIVE CURVATURE (I.E. * TRANS(X)*A*X<0) IF IT EXISTS. * II INF INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. THE * DIRECTION OF NEGATIVE CURVATURE EXISTS ONLY IF INF>0. * * METHOD : * P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND * LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974) * PP. 311-350. * SUBROUTINE MXDPGD(N,A,X,INF) INTEGER N, INF DOUBLE PRECISION A(N*(N+1)/2), X(N) INTEGER I, J, II, IJ DOUBLE PRECISION ZERO, ONE PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) C C RIGHT HAND SIDE FORMATION C DO 1 I=1,N X(I) = ZERO 1 CONTINUE IF (INF .LE. 0) RETURN X(INF) = ONE C C BACK SUBSTITUTION C II = INF*(INF-1)/2 DO 3 I = INF-1, 1, -1 IJ = II DO 2 J = I+1, INF IJ = IJ + J - 1 X(I) = X(I) - A(IJ)*X(J) 2 CONTINUE II = II - I 3 CONTINUE RETURN END * SUBROUTINE MXDPGF ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * FACTORIZATION A+E=L*D*TRANS(L) OF A DENSE SYMMETRIC POSITIVE DEFINITE * MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE DEFINITE MATRICES AND * L IS A LOWER TRIANGULAR MATRIX. IF A IS SUFFICIENTLY POSITIVE * DEFINITE THEN E=0. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) ON INPUT A GIVEN DENSE SYMMETRIC (USUALLY POSITIVE * DEFINITE) MATRIX A STORED IN THE PACKED FORM. ON OUTPUT THE * COMPUTED FACTORIZATION A+E=L*D*TRANS(L). * IO INF AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF * INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF * INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF * INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE * MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION * PROCESS. * RU ALF ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON * OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE * FACTORIZATION PROCESS (IF INF>0). * RO TAU MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E. * * METHOD : * P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND * LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974) * PP. 311-350. * SUBROUTINE MXDPGF(N,A,INF,ALF,TAU) DOUBLE PRECISION ALF,TAU INTEGER INF,N DOUBLE PRECISION A(*) DOUBLE PRECISION BET,DEL,GAM,RHO,SIG,TOL INTEGER I,IJ,IK,J,K,KJ,KK,L L = 0 INF = 0 TOL = ALF * * ESTIMATION OF THE MATRIX NORM * ALF = 0.0D0 BET = 0.0D0 GAM = 0.0D0 TAU = 0.0D0 KK = 0 DO 20 K = 1,N KK = KK + K BET = MAX(BET,ABS(A(KK))) KJ = KK DO 10 J = K + 1,N KJ = KJ + J - 1 GAM = MAX(GAM,ABS(A(KJ))) 10 CONTINUE 20 CONTINUE BET = MAX(TOL,BET,GAM/N) * DEL = TOL*BET DEL = TOL*MAX(BET,1.0D0) KK = 0 DO 60 K = 1,N KK = KK + K * * DETERMINATION OF A DIAGONAL CORRECTION * SIG = A(KK) IF (ALF.GT.SIG) THEN ALF = SIG L = K END IF GAM = 0.0D0 KJ = KK DO 30 J = K + 1,N KJ = KJ + J - 1 GAM = MAX(GAM,ABS(A(KJ))) 30 CONTINUE GAM = GAM*GAM RHO = MAX(ABS(SIG),GAM/BET,DEL) IF (TAU.LT.RHO-SIG) THEN TAU = RHO - SIG INF = -1 END IF * * GAUSSIAN ELIMINATION * A(KK) = RHO KJ = KK DO 50 J = K + 1,N KJ = KJ + J - 1 GAM = A(KJ) A(KJ) = GAM/RHO IK = KK IJ = KJ DO 40 I = K + 1,J IK = IK + I - 1 IJ = IJ + 1 A(IJ) = A(IJ) - A(IK)*GAM 40 CONTINUE 50 CONTINUE 60 CONTINUE IF (L.GT.0 .AND. ABS(ALF).GT.DEL) INF = L RETURN END * SUBROUTINE MXDPGN ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * ESTIMATION OF THE MINIMUM EIGENVALUE AND THE CORRESPONDING EIGENVECTOR * OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE * FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RO X(N) ESTIMATED EIGENVECTOR. * RO ALF ESTIMATED EIGENVALUE. * II JOB OPTION. IF JOB=0 THEN ONLY ESTIMATED EIGENVALUE IS * COMPUTED. IS JOB>0 THEN BOTH ESTIMATED EIGENVALUE AND * ESTIMATED EIGENVECTOR ARE COMPUTED BY JOB ITERATIONS. * * SUBPROGRAMS USED : * S MXDPGB BACK SUBSTITUTION. * RF MXVDOT DOT PRODUCT OF TWO VECTORS. * * METHOD : * A.K.CLINE, C.B.MOLER, G.W.STEWART, J.H.WILKINSON : AN ESTIMATE * FOR THE CONDITION NUMBER OF A MATRIX. SIAM J. NUMER. ANAL. 16 * (1979) 368-373. * SUBROUTINE MXDPGN(N,A,X,ALF,JOB) INTEGER N,JOB DOUBLE PRECISION A(N*(N+1)/2),X(N),ALF DOUBLE PRECISION XP,XM,FP,FM,MXVDOT INTEGER I,K,IK,KK DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) C C COMPUTATION OF THE VECTOR V WITH POSSIBLE MAXIMUM NORM SUCH C THAT L*D**(1/2)*V=U WHERE U HAS ELEMENTS +1 OR -1 C DO 2 I=1,N X(I)=ZERO 2 CONTINUE KK=0 DO 6 K=1,N KK=KK+K XP=-X(K)+ONE XM=-X(K)-ONE FP=ABS(XP) FM=ABS(XM) IK=KK DO 3 I=K+1,N IK=IK+I-1 FP=FP+ABS(X(I)+A(IK)*XP) FM=FM+ABS(X(I)+A(IK)*XM) 3 CONTINUE IF(FP.GE.FM) THEN X(K)=XP IK=KK DO 4 I=K+1,N IK=IK+I-1 X(I)=X(I)+A(IK)*XP 4 CONTINUE ELSE X(K)=XM IK=KK DO 5 I=K+1,N IK=IK+I-1 X(I)=X(I)+A(IK)*XM 5 CONTINUE ENDIF 6 CONTINUE C C COMPUTATION OF THE VECTOR X SUCH THAT C D**(1/2)*TRANS(L)*X=V C FM=ZERO KK=0 DO 7 K=1,N KK=KK+K IF (JOB.LE.0) THEN FP=SQRT(A(KK)) X(K)=X(K)/FP FM=FM+X(K)*X(K) ELSE X(K)=X(K)/A(KK) ENDIF 7 CONTINUE FP=DBLE(N) IF (JOB.LE.0) THEN C C FIRST ESTIMATION OF THE MINIMUM EIGENVALUE BY THE FORMULA C ALF=(TRANS(U)*U)/(TRANS(V)*V) C ALF=FP/FM RETURN ENDIF FM=ZERO KK=N*(N+1)/2 DO 9 K=N,1,-1 IK=KK DO 8 I=K+1,N IK=IK+I-1 X(K)=X(K)-A(IK)*X(I) 8 CONTINUE FM=FM+X(K)*X(K) KK=KK-K 9 CONTINUE FM=SQRT(FM) IF (JOB.LE.1) THEN C C SECOND ESTIMATION OF THE MINIMUM EIGENVALUE BY THE FORMULA C ALF=SQRT(TRANS(U)*U)/SQRT(TRANS(X)*X) C ALF=SQRT(FP)/FM ELSE C C INVERSE ITERATIONS C DO 11 K=2,JOB C C SCALING THE VECTOR X BY ITS NORM C DO 10 I=1,N X(I)=X(I)/FM 10 CONTINUE CALL MXDPGB(N,A,X,0) FM=SQRT(MXVDOT(N,X,X)) 11 CONTINUE ALF=ONE/FM ENDIF C C SCALING THE VECTOR X BY ITS NORM C DO 12 I=1,N X(I)=X(I)/FM 12 CONTINUE RETURN END * FUNCTION MXDPGP ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COMPUTATION OF THE NUMBER MXDPGP=TRANS(X)*D**(-1)*Y WHERE D IS A * DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR MXDPGP COMPUTED NUMBER MXDPGP=TRANS(X)*D**(-1)*Y. * FUNCTION MXDPGP(N,A,X,Y) INTEGER N DOUBLE PRECISION A(*),X(*),Y(*),MXDPGP DOUBLE PRECISION TEMP INTEGER I,J J = 0 TEMP = 0.0D0 DO 10 I = 1,N J = J + I TEMP = TEMP + X(I)*Y(I)/A(J) 10 CONTINUE MXDPGP = TEMP RETURN END * SUBROUTINE MXDPGS ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE * FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RI ALF SCALING FACTOR. * SUBROUTINE MXDPGS(N,A,ALF) DOUBLE PRECISION ALF INTEGER N DOUBLE PRECISION A(*) INTEGER I,J J = 0 DO 10 I = 1,N J = J + I A(J) = A(J)*ALF 10 CONTINUE RETURN END * SUBROUTINE MXDPGU ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E IN THE * FACTORED FORM A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF. * THE CORRECTION IS DEFINED AS A+E:=A+E+ALF*X*TRANS(X) WHERE ALF IS A * GIVEN SCALING FACTOR AND X IS A GIVEN VECTOR. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE MXDPGF. * RI ALF SCALING FACTOR IN THE CORRECTION TERM. * RI X(N) VECTOR IN THE CORRECTION TERM. * RA Y(N) AUXILIARY VECTOR. * * METHOD : * P.E.GILL, W.MURRAY, M.SAUNDERS: METHODS FOR COMPUTING AND MODIFYING * THE LDV FACTORS OF A MATRIX, MATH. OF COMP. 29 (1974) PP. 1051-1077. * SUBROUTINE MXDPGU(N,A,ALF,X,Y) DOUBLE PRECISION ZERO,ONE,FOUR,CON PARAMETER (ZERO=0.0D0,ONE=1.0D0,FOUR=4.0D0,CON=1.0D-8) DOUBLE PRECISION ALF,ALFR INTEGER N DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION B,D,P,R,T,TO INTEGER I,II,IJ,J IF (ALF.GE.ZERO) THEN * * FORWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NONNEGATIVE * ALFR = SQRT(ALF) CALL MXVSCL(N,ALFR,X,Y) TO = ONE II = 0 DO 30 I = 1,N II = II + I D = A(II) P = Y(I) T = TO + P*P/D R = TO/T A(II) = D/R B = P/ (D*T) IF (A(II).LE.FOUR*D) THEN * * AN EASY FORMULA FOR LIMITED DIAGONAL ELEMENT * IJ = II DO 10 J = I + 1,N IJ = IJ + J - 1 D = A(IJ) Y(J) = Y(J) - P*D A(IJ) = D + B*Y(J) 10 CONTINUE ELSE * * A MORE COMPLICATE BUT NUMERICALLY STABLE FORMULA FOR UNLIMITED * DIAGONAL ELEMENT * IJ = II DO 20 J = I + 1,N IJ = IJ + J - 1 D = A(IJ) A(IJ) = R*D + B*Y(J) Y(J) = Y(J) - P*D 20 CONTINUE END IF TO = T 30 CONTINUE ELSE * * BACKWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NEGATIVE * ALFR = SQRT(-ALF) CALL MXVSCL(N,ALFR,X,Y) TO = ONE IJ = 0 DO 50 I = 1,N D = Y(I) DO 40 J = 1,I - 1 IJ = IJ + 1 D = D - A(IJ)*Y(J) 40 CONTINUE Y(I) = D IJ = IJ + 1 TO = TO - D*D/A(IJ) 50 CONTINUE IF (TO.LE.ZERO) TO = CON II = N* (N+1)/2 DO 70 I = N,1,-1 D = A(II) P = Y(I) T = TO + P*P/D A(II) = D*TO/T B = -P/ (D*TO) TO = T IJ = II DO 60 J = I + 1,N IJ = IJ + J - 1 D = A(IJ) A(IJ) = D + B*Y(J) Y(J) = Y(J) + P*D 60 CONTINUE II = II - I 70 CONTINUE END IF RETURN END * SUBROUTINE MXDPRB ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC * POSITIVE DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R. * RU X(N) ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR * EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR * EQUATIONS. * II JOB OPTION. IF JOB=0 THEN X:=A**(-1)*X. IF JOB>0 THEN * X:=TRANS(R)**(-1)*X. IF JOB<0 THEN X:=R**(-1)*X. * * METHOD : * BACK SUBSTITUTION * SUBROUTINE MXDPRB(N,A,X,JOB) INTEGER JOB,N DOUBLE PRECISION A(*),X(*) INTEGER I,II,IJ,J IF (JOB.GE.0) THEN * * PHASE 1 : X:=TRANS(R)**(-1)*X * IJ = 0 DO 20 I = 1,N DO 10 J = 1,I - 1 IJ = IJ + 1 X(I) = X(I) - A(IJ)*X(J) 10 CONTINUE IJ = IJ + 1 X(I) = X(I)/A(IJ) 20 CONTINUE END IF IF (JOB.LE.0) THEN * * PHASE 2 : X:=R**(-1)*X * II = N* (N+1)/2 DO 40 I = N,1,-1 IJ = II DO 30 J = I + 1,N IJ = IJ + J - 1 X(I) = X(I) - A(IJ)*X(J) 30 CONTINUE X(I) = X(I)/A(II) II = II - I 40 CONTINUE END IF RETURN END * SUBROUTINE MXDPRC ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CORRECTION OF A SINGULAR DENSE SYMMETRIC POSITIVE SEMIDEFINITE MATRIX * A DECOMPOSED AS A=TRANS(R)*R. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * IO INF AN INFORMATION OBTAINED IN THE CORRECTION PROCESS. IF * INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE. IF * INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE. * PROCESS. * RI TOL DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. * SUBROUTINE MXDPRC(N,A,INF,TOL) INTEGER N,INF DOUBLE PRECISION A(N*(N+1)/2),TOL DOUBLE PRECISION TOL1,TEMP INTEGER L,I DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) INF=0 TOL1=SQRT(TOL) TEMP=TOL1 DO 3 I=1,N*(N+1)/2 TEMP=MAX(TEMP,ABS(A(I))) 3 CONTINUE TEMP=TEMP*TOL1 L=0 DO 1 I=1,N L=L+I IF (ABS(A(L)).LE.TEMP) THEN A(L)=SIGN(TEMP,A(L)) INF=-1 ENDIF 1 CONTINUE RETURN END * SUBROUTINE MXDPRM ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A GIVEN VECTOR X BY A DENSE SYMMETRIC POSITIVE * DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R. * RU X(N) ON INPUT THE GIVEN VECTOR. ON OUTPUT THE RESULT OF * MULTIPLICATION. * II JOB OPTION. IF JOB=0 THEN X:=A*X. IF JOB>0 THEN X:=R*X. * IF JOB<0 THEN X:=TRANS(R)*X. * SUBROUTINE MXDPRM( N, A, X, JOB) INTEGER N, JOB DOUBLE PRECISION A(N*(N+1)/2), X(N) INTEGER I, J, II, IJ IF (JOB .GE. 0) THEN C C PHASE 1 : X:=R*X C II = 0 DO 3 I = 1, N II = II + I X(I) = A(II) * X(I) IJ = II DO 2 J = I+1, N IJ = IJ + J - 1 X(I) = X(I) + A(IJ) * X(J) 2 CONTINUE 3 CONTINUE ENDIF IF (JOB .LE. 0) THEN C C PHASE 2 : X:=TRANS(R)*X C IJ = N*(N+1)/2 DO 6 I = N, 1, -1 X(I) = A(IJ) * X(I) DO 5 J = I-1, 1, -1 IJ = IJ - 1 X(I) = X(I) + A(IJ) * X(J) 5 CONTINUE IJ = IJ - 1 6 CONTINUE ENDIF RETURN END * SUBROUTINE MXDRGR ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PLANE ROTATION IS APPLIED TO A ROWWISE STORED DENSE RECTANGULAR * MATRIX A. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * RU A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * II K FIRST INDEX OF THE PLANE ROTATION. * II L SECOND INDEX OF THE PLANE ROTATION. * RI CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RI CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * II IER TYPE OF THE PLANE ROTATION. IER=0-GENERAL PLANE ROTATION. * IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED. * * SUBPROGRAMS USED : * S MXVROT PLANE ROTATION APPLIED TO TWO ELEMENTS. * SUBROUTINE MXDRGR(N,A,K,L,CK,CL,IER) INTEGER N,K,L,IER DOUBLE PRECISION A(*),CK,CL INTEGER I,IK,IL IF (IER.NE.0.AND.IER.NE.1) RETURN IK=(K-1)*N IL=(L-1)*N DO 1 I=1,N IK=IK+1 IL=IL+1 CALL MXVROT(A(IK),A(IL),CK,CL,IER) 1 CONTINUE RETURN END * SUBROUTINE MXDRMD ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY * A VECTOR X AND ADDITION OF A SCALED VECTOR ALF*Y. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RI X(N) INPUT VECTOR. * RI ALF SCALING FACTOR. * RI Y(M) INPUT VECTOR. * RO Z(M) OUTPUT VECTOR EQUAL TO A*X+ALF*Y. * SUBROUTINE MXDRMD(N,M,A,X,ALF,Y,Z) INTEGER N,M DOUBLE PRECISION A(M*N),X(N),ALF,Y(M),Z(M) DOUBLE PRECISION TEMP INTEGER I,J,K K=0 DO 2 J=1,M TEMP=ALF*Y(J) DO 1 I=1,N TEMP=TEMP+A(K+I)*X(I) 1 CONTINUE Z(J)=TEMP K=K+N 2 CONTINUE RETURN END * SUBROUTINE MXDRMI ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * ROWWISE STORED DENSE RECTANGULAR MATRIX A IS SET TO BE A PART OF THE * UNIT MATRIX. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RO A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE ONE-DIMENSIONAL * ARRAY. THIS MATRIX IS SET TO TRANS([I,0]). * SUBROUTINE MXDRMI(N,M,A) INTEGER N,M DOUBLE PRECISION A(M*N) INTEGER I,J,K DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) K=0 DO 2 J=1,M DO 1 I=1,N A(I+K)=ZERO IF (I.EQ.J) A(I+K)=ONE 1 CONTINUE K=K+N 2 CONTINUE RETURN END * SUBROUTINE MXDRMM ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY * A VECTOR X. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RI X(N) INPUT VECTOR. * RO Y(M) OUTPUT VECTOR EQUAL TO A*X. * SUBROUTINE MXDRMM(N,M,A,X,Y) INTEGER N,M DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I,J,K K=0 DO 2 J=1,M TEMP=0.0D 0 DO 1 I=1,N TEMP=TEMP+A(K+I)*X(I) 1 CONTINUE Y(J)=TEMP K=K+N 2 CONTINUE RETURN END * FUNCTION MXDRMN ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EUCLIDEAN NORM OF A PART OF THE I-TH COLUMN OF A ROWWISE STORED DENSE * RECTANGULAR MATRIX A IS COMPUTED. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * II I INDEX OF THE COLUMN WHOSE NORM IS COMPUTED. * II J INDEX OF THE FIRST ELEMENT FROM WHICH THE NORM IS COMPUTED. * FUNCTION MXDRMN(N,M,A,I,J) INTEGER N,M,I,J DOUBLE PRECISION A(M*N),MXDRMN DOUBLE PRECISION POM,DEN INTEGER K,L DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) DEN=ZERO L=(J-1)*N DO 1 K=J,M DEN=MAX(DEN,ABS(A(L+I))) L=L+N 1 CONTINUE POM=ZERO IF (DEN.GT.ZERO) THEN L=(J-1)*N DO 2 K=J,M POM=POM+(A(L+I)/DEN)**2 L=L+N 2 CONTINUE ENDIF MXDRMN=DEN*SQRT(POM) RETURN END * SUBROUTINE MXDRMV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * K-TH COLUMN OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A IS COPIED * TO THE VECTOR X. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX A. * II M NUMBER OF ROWS OF THE MATRIX A. * RI A(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RO X(M) OUTPUT VECTOR SUCH THAT X(J)=A(J,K) FOR ALL J. * II K INDEX OF THE ROW BEING COPIED TO THE OUTPUT VECTOR. * SUBROUTINE MXDRMV(N,M,A,X,K) INTEGER N,M,K DOUBLE PRECISION A(*),X(*) INTEGER I,J IF (K.LT.1.OR.K.GT.N) RETURN I=K DO 1 J=1,M X(J)=A(I) I=I+N 1 CONTINUE RETURN END * SUBROUTINE MXDRQF ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * QR DECOMPOSITION OF ROWWISE STORED DENSE RECTANGULAR MATRIX Q USING * HOUSEHOLDER TRANSFORMATIONS WITHOUT PIVOTING. * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX Q. * II M NUMBER OF ROWS OF THE MATRIX Q. * RU Q(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY. * RO R(N*(N+1)/2) UPPER TRIANGULAR MATRIX STORED IN THE PACKED FORM. * * SUBPROGRAMS USED : * S MXDRMN EUCLIDEAN NORM OF A PART OF THE ROWWISE STORED * RECTANGULAR MATRIX COLUMN. * * METHOD : * P.A.BUSSINGER, G.H.GOLUB : LINEAR LEAST SQUARES SOLUTION BY * HOUSEHOLDER TRANSFORMATION. NUMER. MATH. 7 (1965) 269-276. * SUBROUTINE MXDRQF(N,M,Q,R) INTEGER N,M DOUBLE PRECISION Q(M*N),R(N*(N+1)/2) DOUBLE PRECISION ALF,POM,MXDRMN INTEGER I,J,K,L,JP,KP,NM DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) NM=MIN(N,M) C C QR DECOMPOSITION C L=0 KP=0 DO 6 K=1,NM POM=MXDRMN(N,M,Q,K,K) IF (Q(KP+K).NE.ZERO) POM=SIGN(POM,Q(KP+K)) R(L+K)=-POM JP=0 DO 1 J=1,K-1 R(L+J)=Q(JP+K) Q(JP+K)=ZERO JP=JP+N 1 CONTINUE IF (POM.NE.ZERO) THEN C C HOUSEHOLDER TRANSFORMATION C DO 2 J=K,M Q(JP+K)=Q(JP+K)/POM JP=JP+N 2 CONTINUE Q(KP+K)=Q(KP+K)+ONE DO 5 I=K+1,N ALF=ZERO JP=KP DO 3 J=K,M ALF=ALF+Q(JP+K)*Q(JP+I) JP=JP+N 3 CONTINUE ALF=ALF/Q(KP+K) JP=KP DO 4 J=K,M Q(JP+I)=Q(JP+I)-ALF*Q(JP+K) JP=JP+N 4 CONTINUE 5 CONTINUE ENDIF L=L+K KP=KP+N 6 CONTINUE C C EXPLICIT FORMULATION OF THE ORTHOGONAL MATRIX C KP=N*N DO 11 K=N,1,-1 KP=KP-N IF (Q(KP+K).NE.ZERO) THEN DO 9 I=K+1,N ALF=ZERO JP=KP DO 7 J=K,M ALF=ALF+Q(JP+K)*Q(JP+I) JP=JP+N 7 CONTINUE ALF=ALF/Q(KP+K) JP=KP DO 8 J=K,M Q(JP+I)=Q(JP+I)-ALF*Q(JP+K) JP=JP+N 8 CONTINUE 9 CONTINUE JP=KP DO 10 J=K,M Q(JP+K)=-Q(JP+K) JP=JP+N 10 CONTINUE ENDIF Q(KP+K)=Q(KP+K)+ONE 11 CONTINUE RETURN END * SUBROUTINE MXDRQU ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * UPDATE OF A QR DECOMPOSITION. THIS QR DECOMPOSITION IS UPDATED * BY THE RULE Q*R:=Q*R+ALF*X*TRANS(Y). * * PARAMETERS : * II N NUMBER OF COLUMNS OF THE MATRIX Q. * II M NUMBER OF ROWS OF THE MATRIX Q. * RU Q(M*N) RECTANGULAR MATRIX STORED ROWWISE IN THE * ONE-DIMENSIONAL ARRAY (PART OF THE ORTHOGONAL MATRIX). * RU R(N*(N+1)/2) UPPER TRIANGULAR MATRIX STORED IN A PACKED FORM. * RI ALF SCALAR PARAMETER. * RI X(M) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RA Z(N) AUXILIARY VECTOR. * IO INF INFORMATION. IF INF=0 THEN X LIES IN THE COLUMN SPACE OF Q. * IF INF=1 THEN X DOES NOT LIE IN THE COLUMN SPACE OF Q. * * SUBPROGRAMS USED : * RF MXVNOR EUCLIDEAN NORM OF A VECTOR. * S MXVORT COMPUTATION OF THE PLANE ROTATION MATRIX. * S MXVROT PLANE ROTATION IS APPLIED TO TWO NUMBERS. * * METHOD : * J.W.DANIEL, W.B.GRAGG, L.KAUFMAN, G.W.STEWARD : REORTHOGONALIZATION * AND STABLE ALGORITHMS FOR UPDATING THE GRAM-SCHMIDT QR FACTORIZATION. * MATHEMATICS OF COMPUTATION 30 (1976) 772-795. SUBROUTINE MXDRQU(N,M,Q,R,ALF,X,Y,Z,INF) INTEGER N,M,INF DOUBLE PRECISION Q(M*N),R(N*(N+1)/2),ALF,X(M),Y(N),Z(N) DOUBLE PRECISION CK,CL,ZK,ZL,MXVNOR INTEGER J,K,L,KJ,KK,IER DOUBLE PRECISION ONE,CON PARAMETER (ONE=1.0D 0,CON=1.0D-10) INF=0 KK=N*(N+1)/2 C C COMPUTATION OF THE VECTOR TRANS(Q)*X C CALL MXDCMM(N,M,Q,X,Z) IF (M.GT.N) THEN C C IF X DOES NOT LIE IN THE COLUMN SPACE OF Q WE HAVE TO USE C A SUBPROBLEM WHOSE DIMENSION IS BY ONE GREATER (INF=1). C ZK=MXVNOR(M,X) CALL MXDRMD(N,M,Q,Z,-ONE,X,X) ZL=MXVNOR(M,X) IF (ZL.GT.CON*ZK) THEN INF=1 CALL MXVSCL(M,-ONE/ZL,X,X) CALL MXVORT(Z(N),ZL,CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN CALL MXVROT(R(KK),ZL,CK,CL,IER) KJ=N DO 1 J=1,M CALL MXVROT(Q(KJ),X(J),CK,CL,IER) KJ=KJ+N 1 CONTINUE ENDIF ENDIF ENDIF C C APPLICATION OF PLANE ROTATIONS TO THE VECTOR Z SO THAT C TRANS(Q1)*Z=E1 WHERE Q1 IS AN ORTHOGONAL MATRIX (ACCUMULATION OF C THE PLANE ROTATIONS) AND E1 IS THE FIRST COLUMN OF THE UNIT C MATRIX. AT THE SAME TIME BOTH THE UPPER HESSENBERG MATRIX C TRANS(Q1)*R AND THE ORTHOGONAL MATRIX Q*Q1 ARE CONSTRUCTED SO THAT C Q*Q1*R1=Q*Q1*(TRANS(Q1)*R+ALF*E1*TRANS(Y)) WHERE R1 IS AN UPPER C HESSENBERG MATRIX. C DO 4 L=N,2,-1 K=L-1 KK=KK-L CALL MXVORT(Z(K),Z(L),CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN CALL MXVROT(R(KK),Z(L),CK,CL,IER) KJ=KK DO 2 J=L,N KJ=KJ+J-1 CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER) 2 CONTINUE KJ=K DO 3 J=1,M CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER) KJ=KJ+N 3 CONTINUE ENDIF 4 CONTINUE Z(1)=ALF*Z(1) KJ=1 DO 5 J=1,N R(KJ)=R(KJ)+Z(1)*Y(J) KJ=KJ+J 5 CONTINUE C C APPLICATION OF PLANE ROTATIONS TO THE UPPER HESSENBERG MATRIX R1 C GIVEN ABOVE SO THAT R2=TRANS(Q2)*R1 WHERE Q2 IS AN ORTHOGONAL C MATRIX (ACCUMULATION OF THE PLANE ROTATIONS) AND R2 IS AN UPPER C TRIANGULAR MATRIX. WE OBTAIN THE NEW QR DECOMPOSITION Q*Q1*Q2*R2. C KK=1 DO 8 L=2,N K=L-1 CALL MXVORT(R(KK),Z(L),CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN KJ=KK DO 6 J=L,N KJ=KJ+J-1 CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER) 6 CONTINUE KJ=K DO 7 J=1,M CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER) KJ=KJ+N 7 CONTINUE ENDIF KK=KK+L 8 CONTINUE C C BACK TRANSFORMATION OF THE GREATER SUBPROBLEM IF INF=1. C IF (INF.EQ.1) THEN CALL MXVORT(R(KK),ZL,CK,CL,IER) IF (IER.EQ.0.OR.IER.EQ.1) THEN KJ=N DO 9 J=1,M CALL MXVROT(Q(KJ),X(J),CK,CL,IER) KJ=KJ+N 9 CONTINUE ENDIF ENDIF RETURN END * SUBROUTINE MXDSDA ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * A DENSE SYMMETRIC MATRIX A IS AUGMENTED BY THE SCALED UNIT MATRIX * SUCH THAT A:=A+ALF*I (I IS THE UNIT MATRIX OF ORDER N). * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RI ALF SCALING FACTOR. * SUBROUTINE MXDSDA(N,A,ALF) INTEGER N DOUBLE PRECISION A(*),ALF INTEGER I,J J = 0 DO 1 I = 1, N J = J + I A(J)=A(J)+ALF 1 CONTINUE RETURN END * FUNCTION MXDSDL ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF THE MINIMUM DIAGONAL ELEMENT OF A DENSE SYMMETRIC * MATRIX. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * IO INF INDEX OF THE MINIMUM DIAGONAL ELEMENT OF THE MATRIX A. * RR MXDSDL MINIMUM DIAGONAL ELEMENT OF THE MATRIX A. * FUNCTION MXDSDL(N,A,INF) INTEGER N,INF DOUBLE PRECISION A(N*(N+1)/2),MXDSDL DOUBLE PRECISION TEMP INTEGER I,J J=1 INF=1 TEMP=A(1) DO 1 I=2,N J=J+I IF (TEMP.GT.A(J)) THEN INF=J TEMP=A(J) ENDIF 1 CONTINUE MXDSDL=TEMP RETURN END * SUBROUTINE MXDSMA ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DENSE SYMMETRIC MATRIX AUGMENTED BY THE SCALED DENSE SYMMETRIC MATRIX. * * PARAMETERS : * II N ORDER OF THE MATRICES. * RI ALF SCALING FACTOR. * RI A(N*(N+1)/2) INPUT MATRIX. * RI B(N*(N+1)/2) INPUT MATRIX. * RO C(N*(N+1)/2) OUTPUT MATRIX WHERE C:=B+ALF*A. * SUBROUTINE MXDSMA(N,ALF,A,B,C) INTEGER N DOUBLE PRECISION ALF,A(N*(N+1)/2),B(N*(N+1)/2),C(N*(N+1)/2) INTEGER I DO 1 I=1,N*(N+1)/2 C(I)=B(I)+ALF*A(I) 1 CONTINUE RETURN END * SUBROUTINE MXDSMC ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COPYING OF A DENSE SYMMETRIC MATRIX. * * PARAMETERS : * II N ORDER OF THE MATRICES A AND B. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RO B(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM * WHERE B:=A. * SUBROUTINE MXDSMC(N,A,B) INTEGER N DOUBLE PRECISION A(N*(N+1)/2), B(N*(N+1)/2) INTEGER M,I M = N * (N+1)/2 DO 1 I = 1, M B(I) = A(I) 1 CONTINUE RETURN END * SUBROUTINE MXDSMG ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * GERSHGORIN BOUNDS FOR EIGENVALUES OF A DENSE SYMMETRIC MATRIX. * AMIN .LE. ANY EIGENVALUE OF A .LE. AMAX. * * PARAMETERS : * II N DIMENSION OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RO AMIN LOWER BOUND FOR EIGENVALUES OF A. * RO AMAX UPPER BOUND FOR EIGENVALUES OF A. * SUBROUTINE MXDSMG(N,A,AMIN,AMAX) INTEGER N DOUBLE PRECISION A(N*(N+1)/2),AMIN,AMAX DOUBLE PRECISION TEMP INTEGER I,J,K,L DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) AMAX=A(1) AMIN=A(1) K=0 DO 3 I=1,N TEMP=ZERO L=K DO 1 J=1,I-1 L=L+1 TEMP=TEMP+ABS(A(L)) 1 CONTINUE L=L+1 DO 2 J=I+1,N L=L+J-1 TEMP=TEMP+ABS(A(L)) 2 CONTINUE K=K+I AMAX=MAX(AMAX,A(K)+TEMP) AMIN=MIN(AMIN,A(K)-TEMP) 3 CONTINUE RETURN END * SUBROUTINE MXDSMI ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DENSE SYMMETRIC MATRIX A IS SET TO THE UNIT MATRIX WITH THE SAME * ORDER. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RO A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM * WHICH IS SET TO THE UNIT MATRIX (I.E. A:=I). * SUBROUTINE MXDSMI(N,A) INTEGER N DOUBLE PRECISION A(*) INTEGER I,M M = N* (N+1)/2 DO 10 I = 1,M A(I) = 0.0D0 10 CONTINUE M = 0 DO 20 I = 1,N M = M + I A(M) = 1.0D0 20 CONTINUE RETURN END * SUBROUTINE MXDSMM ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * MULTIPLICATION OF A DENSE SYMMETRIC MATRIX A BY A VECTOR X. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR EQUAL TO A*X. * SUBROUTINE MXDSMM(N,A,X,Y) INTEGER N DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I,J,K,L K = 0 DO 30 I = 1,N TEMP = 0.0D0 L = K DO 10 J = 1,I L = L + 1 TEMP = TEMP + A(L)*X(J) 10 CONTINUE DO 20 J = I + 1,N L = L + J - 1 TEMP = TEMP + A(L)*X(J) 20 CONTINUE Y(I) = TEMP K = K + I 30 CONTINUE RETURN END * FUNCTION MXDSMQ ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VALUE OF A QUADRATIC FORM WITH A DENSE SYMMETRIC MATRIX A. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RI X(N) GIVEN VECTOR. * RI Y(N) GIVEN VECTOR. * RR MXDSMQ VALUE OF THE QUADRATIC FORM MXDSMQ=TRANS(X)*A*Y. * DOUBLE PRECISION + FUNCTION MXDSMQ(N,A,X,Y) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER N DOUBLE PRECISION A(*),X(*),Y(*) DOUBLE PRECISION TEMP,TEMP1,TEMP2 INTEGER I,J,K TEMP = ZERO K = 0 DO 20 I = 1,N TEMP1 = ZERO TEMP2 = ZERO DO 10 J = 1,I - 1 K = K + 1 TEMP1 = TEMP1 + A(K)*X(J) TEMP2 = TEMP2 + A(K)*Y(J) 10 CONTINUE K = K + 1 TEMP = TEMP + X(I)* (TEMP2+A(K)*Y(I)) + Y(I)*TEMP1 20 CONTINUE MXDSMQ = TEMP RETURN END * SUBROUTINE MXDSMR ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PLANE ROTATION IS APPLIED TO A DENSE SYMMETRIC MATRIX A. THE CASE * K=L+1 IS REQUIRED. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * II K FIRST INDEX OF PLANE ROTATION. * II L SECOND INDEX OF PLANE ROTATION. * RO CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RO CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * IO IER INFORMATION ON THE TRANSFORMATION. IER<0-K OR L OUT OF * RANGE. IER=0-PLANE ROTATION. IER=1-PERMUTATION. * IER=2-TRANSFORMATION SUPPRESSED. * * SUBPROGRAMS USED : * S MXVROT PLANE ROTATION IS APPLIED TO TWO NUMBERS. * SUBROUTINE MXDSMR(N,A,K,L,CK,CL,IER) INTEGER N,K,L,IER DOUBLE PRECISION A(*),CK,CL DOUBLE PRECISION AKK,AKL,ALL,CKK,CKL,CLL INTEGER J,KJ,LJ,KK,KL,LL IF (IER.NE.0.AND.IER.NE.1) RETURN IF (K.NE.L+1) THEN IER=-1 RETURN ENDIF LJ=L*(L-1)/2 DO 1 J=1,N IF (J.LE.L) THEN LJ=LJ+1 KJ=LJ+L ELSE LJ=LJ+J-1 KJ=LJ+1 ENDIF IF (J.NE.K.AND.J.NE.L) THEN CALL MXVROT(A(KJ),A(LJ),CK,CL,IER) ENDIF 1 CONTINUE IF (IER.EQ.0) THEN CKK=CK**2 CKL=CK*CL CLL=CL**2 LL=L*(L+1)/2 KL=LL+L KK=LL+K AKL=(CKL+CKL)*A(KL) AKK=CKK*A(KK)+CLL*A(LL)+AKL ALL=CLL*A(KK)+CKK*A(LL)-AKL A(KL)=(CLL-CKK)*A(KL)+CKL*(A(KK)-A(LL)) A(KK)=AKK A(LL)=ALL ELSE LL=L*(L+1)/2 KK=LL+K AKK=A(KK) A(KK)=A(LL) A(LL)=AKK ENDIF RETURN END * SUBROUTINE MXDSMS ALL SYSTEMS 91/12/01 C PORTABILITY : ALL SYSTEMS C 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALING OF A DENSE SYMMETRIC MATRIX. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM * WHICH IS SCALED BY THE VALUE ALF (I.E. A:=ALF*A). * RI ALF SCALING FACTOR. * SUBROUTINE MXDSMS( N, A, ALF) INTEGER N DOUBLE PRECISION A(N*(N+1)/2), ALF INTEGER I,M M = N * (N+1) / 2 DO 1 I = 1, M A(I) = A(I) * ALF 1 CONTINUE RETURN END * SUBROUTINE MXDSMU ALL SYSTEMS 89/12/01 C PORTABILITY : ALL SYSTEMS C 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * UPDATE OF A DENSE SYMMETRIC MATRIX A. THIS UPDATE IS DEFINED AS * A:=A+ALF*X*TRANS(X) WHERE ALF IS A GIVEN SCALING FACTOR AND X IS * A GIVEN VECTOR. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RU A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RI ALF SCALING FACTOR IN THE CORRECTION TERM. * RI X(N) VECTOR IN THE CORRECTION TERM. * SUBROUTINE MXDSMU( N, A, ALF, X) INTEGER N DOUBLE PRECISION A(N*(N+1)/2), X(N), ALF DOUBLE PRECISION TEMP INTEGER I, J, K K = 0 DO 2 I = 1, N TEMP = ALF * X(I) DO 1 J = 1, I K = K + 1 A(K) = A(K) + TEMP * X(J) 1 CONTINUE 2 CONTINUE RETURN END * SUBROUTINE MXDSMV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * K-TH ROW OF A DENSE SYMMETRIC MATRIX A IS COPIED TO THE VECTOR X. * * PARAMETERS : * II N ORDER OF THE MATRIX A. * RI A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM. * RO X(N) OUTPUT VECTOR. * II K INDEX OF COPIED ROW. * SUBROUTINE MXDSMV(N,A,X,K) INTEGER K,N DOUBLE PRECISION A(*),X(*) INTEGER I,L L = K* (K-1)/2 DO 10 I = 1,N IF (I.LE.K) THEN L = L + 1 ELSE L = L + I - 1 END IF X(I) = A(L) 10 CONTINUE RETURN END * SUBROUTINE MXVCOP ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * COPYING OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:= X. * SUBROUTINE MXVCOP(N,X,Y) INTEGER N DOUBLE PRECISION X(*),Y(*) INTEGER I DO 10 I = 1,N Y(I) = X(I) 10 CONTINUE RETURN END * SUBROUTINE MXVDIF ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTOR DIFFERENCE. * * PARAMETERS : * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= X - Y. * SUBROUTINE MXVDIF(N,X,Y,Z) INTEGER N DOUBLE PRECISION X(*),Y(*),Z(*) INTEGER I DO 10 I = 1,N Z(I) = X(I) - Y(I) 10 CONTINUE RETURN END * SUBROUTINE MXVDIR ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTOR AUGMENTED BY THE SCALED VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI A SCALING FACTOR. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= Y + A*X. * SUBROUTINE MXVDIR(N,A,X,Y,Z) DOUBLE PRECISION A INTEGER N DOUBLE PRECISION X(*),Y(*),Z(*) INTEGER I DO 10 I = 1,N Z(I) = Y(I) + A*X(I) 10 CONTINUE RETURN END * FUNCTION MXVDOT ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DOT PRODUCT OF TWO VECTORS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR MXVDOT VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y. * FUNCTION MXVDOT(N,X,Y) INTEGER N DOUBLE PRECISION X(*),Y(*),MXVDOT DOUBLE PRECISION TEMP INTEGER I TEMP = 0.0D0 DO 10 I = 1,N TEMP = TEMP + X(I)*Y(I) 10 CONTINUE MXVDOT = TEMP RETURN END * SUBROUTINE MXVINA ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES. * * PARAMETERS : * II N DIMENSION OF THE INTEGER VECTOR. * IU IX(N) INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I)) * FOR ALL I. * SUBROUTINE MXVINA(N,IX) INTEGER N INTEGER IX(*) INTEGER I DO 10 I = 1,N IX(I) = ABS(IX(I)) IF (IX(I).GT.10) IX(I) = IX(I) - 10 10 CONTINUE RETURN END * SUBROUTINE MXVIND ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION. * * PARAMETERS : * IU IX(N) INTEGER VECTOR. * II I INDEX OF THE CHANGED ELEMENT. * II JOB CHANGE SPECIFICATION. IS JOB.EQ.0 THEN IX(I)=10-IX(I). * SUBROUTINE MXVIND(IX,I,JOB) INTEGER IX(*),I,JOB IF (JOB.EQ.0) IX(I)=10-IX(I) RETURN END * SUBROUTINE MXVINS ALL SYSTEMS 90/12/01 C PORTABILITY : ALL SYSTEMS C 90/12/01 LU : ORIGINAL VERSION * * PURPOSE : * INITIATION OF THE INTEGER VECTOR. * * PARAMETERS : * II N DIMENSION OF THE INTEGER VECTOR. * II IP INTEGER PARAMETER. * IO IX(N) INTEGER VECTOR SUCH THAT IX(I)=IP FOR ALL I. * SUBROUTINE MXVINS(N,IP,IX) INTEGER N,IP,IX(N) INTEGER I DO 1 I=1,N IX(I)=IP 1 CONTINUE RETURN END * SUBROUTINE MXVINV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION. * * PARAMETERS : * II N VECTOR DIMENSION. * IU IX(N) INTEGER VECTOR. * II I INDEX OF THE CHANGED ELEMENT. * II JOB CHANGE SPECIFICATION * SUBROUTINE MXVINV(IX,I,JOB) INTEGER I,JOB INTEGER IX(*) IF ((IX(I).EQ.3.OR.IX(I).EQ.5) .AND. JOB.LT.0) IX(I) = IX(I) + 1 IF ((IX(I).EQ.4.OR.IX(I).EQ.6) .AND. JOB.GT.0) IX(I) = IX(I) - 1 IX(I) = -IX(I) RETURN END * FUNCTION MXVMAX ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * L-INFINITY NORM OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RR MXVMAX L-INFINITY NORM OF THE VECTOR X. * FUNCTION MXVMAX(N,X) INTEGER N DOUBLE PRECISION X(*),MXVMAX INTEGER I MXVMAX = 0.0D0 DO 10 I = 1,N MXVMAX = MAX(MXVMAX,ABS(X(I))) 10 CONTINUE RETURN END * SUBROUTINE MXVNEG ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * CHANGE THE SIGNS OF VECTOR ELEMENTS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:= - X. * SUBROUTINE MXVNEG(N,X,Y) INTEGER N DOUBLE PRECISION X(*),Y(*) INTEGER I DO 10 I = 1,N Y(I) = -X(I) 10 CONTINUE RETURN END * FUNCTION MXVNOR ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * EUCLIDEAN NORM OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RR MXVNOR EUCLIDEAN NORM OF X. * FUNCTION MXVNOR(N,X) INTEGER N DOUBLE PRECISION X(*),MXVNOR DOUBLE PRECISION POM,DEN INTEGER I DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D 0) DEN=ZERO DO 1 I=1,N DEN=MAX(DEN,ABS(X(I))) 1 CONTINUE POM=ZERO IF (DEN.GT.ZERO) THEN DO 2 I=1,N POM=POM+(X(I)/DEN)**2 2 CONTINUE ENDIF MXVNOR=DEN*SQRT(POM) RETURN END * SUBROUTINE MXVORT ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR PLANE ROTATION. * * PARAMETERS : * RU XK FIRST VALUE FOR PLANE ROTATION (XK IS TRANSFORMED TO * SQRT(XK**2+XL**2)) * RU XL SECOND VALUE FOR PLANE ROTATION (XL IS TRANSFORMED TO * ZERO) * RO CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RO CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * IO IER INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE * ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED. * SUBROUTINE MXVORT(XK,XL,CK,CL,IER) DOUBLE PRECISION CK,CL,XK,XL INTEGER IER DOUBLE PRECISION DEN,POM IF (XL.EQ.0.0D0) THEN IER = 2 ELSE IF (XK.EQ.0.0D0) THEN XK = XL XL = 0.0D0 IER = 1 ELSE IF (ABS(XK).GE.ABS(XL)) THEN POM = XL/XK DEN = SQRT(1.0D0+POM*POM) CK = 1.0D0/DEN CL = POM/DEN XK = XK*DEN ELSE POM = XK/XL DEN = SQRT(1.0D0+POM*POM) CL = 1.0D0/DEN CK = POM/DEN XK = XL*DEN END IF XL = 0.0D0 IER = 0 END IF RETURN END * SUBROUTINE MXVROT ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * PLANE ROTATION IS APPLIED TO TWO VALUES. * * PARAMETERS : * RU XK FIRST VALUE FOR PLANE ROTATION. * RU XL SECOND VALUE FOR PLANE ROTATION. * RI CK DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * RI CL OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX. * II IER INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE * ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED. * SUBROUTINE MXVROT(XK,XL,CK,CL,IER) DOUBLE PRECISION CK,CL,XK,XL INTEGER IER DOUBLE PRECISION YK,YL IF (IER.EQ.0) THEN YK = XK YL = XL XK = CK*YK + CL*YL XL = CL*YK - CK*YL ELSE IF (IER.EQ.1) THEN YK = XK XK = XL XL = YK END IF RETURN END * SUBROUTINE MXVSAV ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RU Y(N) UPDATE VECTOR WHERE Y:= X - Y. * SUBROUTINE MXVSAV(N,X,Y) INTEGER N DOUBLE PRECISION X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I DO 10 I = 1,N TEMP = Y(I) Y(I) = X(I) - Y(I) X(I) = TEMP 10 CONTINUE RETURN END * SUBROUTINE MXVSCL ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SCALING OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI A SCALING FACTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:= A*X. * SUBROUTINE MXVSCL(N,A,X,Y) DOUBLE PRECISION A INTEGER N DOUBLE PRECISION X(*),Y(*) INTEGER I DO 10 I = 1,N Y(I) = A*X(I) 10 CONTINUE RETURN END * SUBROUTINE MXVSET ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI A INITIAL VALUE. * RO X(N) OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I. * SUBROUTINE MXVSET(N,A,X) DOUBLE PRECISION A INTEGER N DOUBLE PRECISION X(*) INTEGER I DO 10 I = 1,N X(I) = A 10 CONTINUE RETURN END * SUBROUTINE MXVSUM ALL SYSTEMS 88/12/01 * PORTABILITY : ALL SYSTEMS * 88/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SUM OF TWO VECTORS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= X + Y. * SUBROUTINE MXVSUM(N,X,Y,Z) INTEGER N DOUBLE PRECISION X(*),Y(*),Z(*) INTEGER I DO 10 I = 1,N Z(I) = X(I) + Y(I) 10 CONTINUE RETURN END