* 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 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 MXDPGI ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * INVERSION OF A DENSE SYMMETRIC MATRIX A+E USING THE DECOMPOSITION * A+E=L*D*TRANS(D) OBTAINED BY THE SUBROUTINE MXDPGF. * * PARAMETERS : * II N ORDER OF THE MATRIX A * RU A(N*(N+1)/2) ON INPUT THE DECOMPOSITION A+E=L*D*TRANS(L) * OBTAINED BY THE SUBROUTINE MXDPGF. * ON OUTPUT THE INVERSION (A+E)**(-1). * * METHOD : * INVERSION OF THE LOWER TRIANGULAR MATRIX L AND BACK MULTIPLICATION * (A+E)**(-1) = TRANS(L)**(-1)*D**(-1)*L**(-1). * SUBROUTINE MXDPGI(N,A) INTEGER N DOUBLE PRECISION A(*) DOUBLE PRECISION AII,AIJ,AIK INTEGER I,J,K,II,IJ,JJ,IK,KJ,KK * * INVERSION OF THE LOWER TRIANGULAR MATRIX L * II=0 DO 3 I=1,N II=II+I A(II)=1.0D 0/A(II) IJ=II DO 2 J=I+1,N IJ=IJ+J-1 AIJ=-A(IJ) IK=II KJ=IJ DO 1 K=I+1,J-1 IK=IK+K-1 KJ=KJ+1 AIJ=AIJ-A(IK)*A(KJ) 1 CONTINUE A(IJ)=AIJ 2 CONTINUE 3 CONTINUE * * BACK MULTIPLICATION (A+E)**(-1)= TRANS(L**(-1))*D**(-1)*L**(-1) * II=0 DO 7 I=1,N II=II+I AII=A(II) IK=II KK=II DO 4 K=I+1,N IK=IK+K-1 KK=KK+K AIK=A(IK)*A(KK) AII=AII+A(IK)*AIK A(IK)=AIK 4 CONTINUE A(II)=AII IJ=II JJ=II DO 6 J=I+1,N IJ=IJ+J-1 JJ=JJ+J AIJ=A(IJ) IK=IJ KJ=JJ DO 5 K=J+1,N IK=IK+K-1 KJ=KJ+K-1 AIJ=AIJ+A(IK)*A(KJ) 5 CONTINUE A(IJ)=AIJ 6 CONTINUE 7 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 UXDPGP=TRANS(X)*D**(-1)*Y WHERE D IS A * DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE * SUBROUTINE UXDPGF. * * 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 UXDPGF. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR MXDPGP COMPUTED NUMBER UXDPGP=TRANS(X)*D**(-1)*Y. * FUNCTION MXDPGP(N,A,X,Y) INTEGER N DOUBLE PRECISION A(N*(N+1)/2), X(N), Y(N), MXDPGP DOUBLE PRECISION TEMP INTEGER I,J J = 0 TEMP = 0.0D0 DO 1 I = 1, N J = J + I TEMP = TEMP + X(I)*Y(I)/A(J) 1 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 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. * II M NUMBER OF ROWS 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 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 * 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 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 * 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 END IF 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 END IF IF (J.NE.K.AND.J.NE.L) THEN CALL MXVROT(A(KJ),A(LJ),CK,CL,IER) END IF 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 END IF RETURN END * SUBROUTINE MXDSMS ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 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(*),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 * PORTABILITY : ALL SYSTEMS * 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(*), X(*), 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 * FUNCTION MXVDEL ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * SQUARED NORM OF A SHIFTED VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI A SCALING FACTOR. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR UXVDEL SQUARED NORM OF Y+A*X. * FUNCTION MXVDEL(N,A,X,Y) INTEGER N DOUBLE PRECISION A, X(*), Y(*), MXVDEL INTEGER I DOUBLE PRECISION TEMP TEMP=0.0D 0 DO 1 I=1,N TEMP=TEMP+(Y(I)+A*X(I))**2 1 CONTINUE MXVDEL=TEMP 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. * DOUBLE PRECISION + FUNCTION MXVDOT(N,X,Y) INTEGER N DOUBLE PRECISION X(*),Y(*) 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 * PORTABILITY : ALL SYSTEMS * 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(*) 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 * SUBROUTINE MXVLIN ALL SYSTEMS 92/12/01 * PORTABILITY : ALL SYSTEMS * 92/12/01 LU : ORIGINAL VERSION * * PURPOSE : * LINEAR COMBINATION OF TWO VECTORS. * * PARAMETERS : * II N VECTOR DIMENSION. * RI A SCALING FACTOR. * RI X(N) INPUT VECTOR. * RI B SCALING FACTOR. * RI Y(N) INPUT VECTOR. * RO Z(N) OUTPUT VECTOR WHERE Z:= A*X + B*Y. * SUBROUTINE MXVLIN(N,A,X,B,Y,Z) INTEGER N DOUBLE PRECISION A, X(*), B, Y(*), Z(*) INTEGER I DO 1 I=1,N Z(I) = A*X(I) + B*Y(I) 1 CONTINUE 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. * DOUBLE PRECISION + FUNCTION MXVMAX(N,X) INTEGER N DOUBLE PRECISION X(*) INTEGER I MXVMAX = 0.0D0 DO 10 I = 1,N MXVMAX = MAX(MXVMAX,ABS(X(I))) 10 CONTINUE RETURN END * FUNCTION MXVMX2 ALL SYSTEMS 95/12/01 * PORTABILITY : ALL SYSTEMS * 95/12/01 SI : ORIGINAL VERSION * * PURPOSE : * L-INFINITY NORM OF VECTOR DIFFERENCE. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR UXVMX2 L-INFINITY NORM OF X-Y. * FUNCTION MXVMX2(N,X,Y) INTEGER N DOUBLE PRECISION X(*),Y(*),MXVMX2 INTEGER I DOUBLE PRECISION TEMP TEMP=0.0D 0 DO 10 I=1,N TEMP=MAX(TEMP,ABS(X(I)-Y(I))) 10 CONTINUE MXVMX2=TEMP END * SUBROUTINE MXVMUL ALL SYSTEMS 89/12/01 * PORTABILITY : ALL SYSTEMS * 89/12/01 LU : ORIGINAL VERSION * * PURPOSE : * VECTOR IS PREMULTIPLIED BY THE POWER OF A DIAGONAL MATRIX. * * PARAMETERS : * II N VECTOR DIMENSION. * RI D(N) DIAGONAL MATRIX STORED AS A VECTOR WITH N ELEMENTS. * RI X(N) INPUT VECTOR. * RO Y(N) OUTPUT VECTOR WHERE Y:=(D**K)*X. * II K INTEGER EXPONENT. * SUBROUTINE MXVMUL(N,D,X,Y,K) INTEGER K,N DOUBLE PRECISION D(*),X(*),Y(*) INTEGER I IF (K.EQ.0) THEN CALL MXVCOP(N,X,Y) ELSE IF (K.EQ.1) THEN DO 10 I = 1,N Y(I) = X(I)*D(I) 10 CONTINUE ELSE IF (K.EQ.-1) THEN DO 20 I = 1,N Y(I) = X(I)/D(I) 20 CONTINUE ELSE DO 30 I = 1,N Y(I) = X(I)*D(I)**K 30 CONTINUE END IF 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 MXVNM2 ALL SYSTEMS 90/12/01 * PORTABILITY : ALL SYSTEMS * 90/12/01 SI : ORIGINAL VERSION * * PURPOSE : * EUCLIDEAN NORM OF VECTOR DIFFERENCE. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RI Y(N) INPUT VECTOR. * RR MXVNM2 EUCLIDEAN NORM OF X-Y. * DOUBLE PRECISION + FUNCTION MXVNM2(N,X,Y) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER N DOUBLE PRECISION X(*),Y(*) DOUBLE PRECISION TEMP INTEGER I TEMP = ZERO DO 10 I = 1,N TEMP = TEMP + (X(I)-Y(I))**2 10 CONTINUE MXVNM2 = SQRT(TEMP) 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 DEN=0.0D 0 DO 1 I=1,N DEN=MAX(DEN,ABS(X(I))) 1 CONTINUE POM=0.0D 0 IF (DEN.GT.0.0D 0) THEN DO 2 I=1,N POM=POM+(X(I)/DEN)**2 2 CONTINUE END IF 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 MXVSAB ALL SYSTEMS 91/12/01 * PORTABILITY : ALL SYSTEMS * 91/12/01 LU : ORIGINAL VERSION * * PURPOSE : * L-1 NORM OF A VECTOR. * * PARAMETERS : * II N VECTOR DIMENSION. * RI X(N) INPUT VECTOR. * RR MXVSAB L-1 NORM OF THE VECTOR X. * FUNCTION MXVSAB(N,X) INTEGER N DOUBLE PRECISION X(*),MXVSAB INTEGER I MXVSAB=0.0D 0 DO 1 I=1,N MXVSAB=MXVSAB+ABS(X(I)) 1 CONTINUE 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