’ÿSGECO FOR.ƒSGEFA FOR1‘ûHILGEN FORH ½ñCODR FORR8åSGESL FORj4SLDR FORƒˆGDIDR FOR (ÓÒSGEDI FORÈ/BLAS FORåPöOUT FOR5¶ÿÁÆo                ¹oÆo SUBROUTINE SGECO(A,LDA,N,IPVT,RCOND,Z) INTEGER LDA,N,IPVT(1) REAL A(LDA,1),Z(1) REAL RCOND C C SGECO FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, SGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW SGECO BY SGESL. C TO COMPUTE INVERSE(A)*C , FOLLOW SGECO BY SGESL. C TO COMPUTE DETERMINANT(A) , FOLLOW SGECO BY SGEDI. C TO COMPUTE INVERSE(A) , FOLLOW SGECO BY SGEDI. C C ON ENTRY C C A REAL(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND REAL C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND .EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z REAL(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C LINPACK SGEFA C BLAS SAXPY,SDOT,SSCAL,SASUM C FORTRAN ABS,AMAX1,SIGN C C INTERNAL VARIABLES C REAL SDOT,EK,T,WK,WKM REAL ANORM,S,SASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C C C COMPUTE 1-NORM OF A C ANORM = 0.0E0 DO 10 J = 1, N ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1)) 10 CONTINUE C C FACTOR C CALL SGEFA(A,LDA,N,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE ?*Z = Y AND TRANS(A)*Y = E . C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID C OVERFLOW. C C SOLVE TRANS(U)*W = E C EK = 1.0E0 DO 20 J = 1, N Z(J) = 0.0E0 20 CONTINUE DO 100 K = 1, N IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,-Z(K)) IF (ABS(EK-Z(K)) .LE. ABS(A(K,K))) GO TO 30 S = ABS(A(K,K))/ABS(EK-Z(K)) CALL SSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = ABS(WK) SM = ABS(WKM) IF (A(K,K) .EQ. 0.0E0) GO TO 40 WK = WK/A(K,K) WKM = WKM/A(K,K) GO TO 50 40 CONTINUE WK = 1.0E0 WKM = 1.0E0 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + ABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + ABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 T = WKM - WK WK = WKM DO 70 J = KP1, N Z(J) = Z(J) + T*A(K,J) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C C SOLVE TRANS(L)*Y = W C DO 120 KB = 1, N K = N + 1 - KB IF (K .LT. N) Z(K) = Z(K) + SDOT(N-K,A(K+1,K),1,Z(K+1),1) IF (ABS(Z(K)) .LE. 1.0E0) GO TO 110 S = 1.0E0/ABS(Z(K)) CALL SSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) C YNORM = 1.0E0 C C SOLVE L*V = Y C DO 140 K = 1, N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K .LT. N) CALL SAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (ABS(Z(K)) .LE. 1.0E0) GO TO 130 S = 1.0E0/ABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = V C DO 160 KB = 1, N K = N + 1 - KB IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 150 S = ABS(A(K,K))/ABS(Z(K)) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K) .NE. 0.0E0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0E0) Z(K) = 1.0E0 T = -Z(K) CALL SAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0E0/SASUM(N,Z,1) CALL SSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0 RETURN END  SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(1),INFO REAL A(LDA,1) C C SGEFA FACTORS A REAL MATRIX BY GAUSSIAN ELIMINATION. C C SGEFA IS USUALLY CALLED BY SGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR SGECO) = (1 + 9/N)*(TIME FOR SGEFA) . C C ON ENTRY C C A REAL(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT SGESL OR SGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN SGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SSCAL,ISAMAX C C INTERNAL VARIABLES C REAL T INTEGER ISAMAX,J,K,KP1,L,NM1 C C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = ISAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (A(L,K) .EQ. 0.0E0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0E0/A(K,K) CALL SSCAL(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL SAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. 0.0E0) INFO = N RETURN END COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730 C C PROGRAM: SYSGEN C C VERSION: 1.0 DATE: 02/22/85 C MICROSOFT FORTRAN C C DESCRIPTION: C C GENERATE A TEST MATRIX AND TEST VECTOR USING HILBERT C COEFFICIENTS AND SET VECTOR TO MATRIX DIAGONAL C ELEMENTS. THE MATRIX IS SYMMETRIC AND POORLY CON- C DITIONED. C C AUTHOR: C C ADAM FRITZ C 133 MAIN STREET C AFTON, NEW YORK 13730 C SUBROUTINE SYSGEN(A, LDA, N, B) C INTEGER LDA, N, I, J, NM1, IP1 REAL A(LDA,1), T, AIJ, B(1) C IF (N .LT. 1 .OR. N .GT. LDA) STOP NM1 = N - 1 IF (N .LT. 2) GO TO 30 DO 20 I=1, NM1 T = 1.E0/FLOAT(2*I-1) A(I,I) = T IP1 = I+1 DO 10 J=IP1, N AIJ = 1.E0/FLOAT(I+J-1) A(I,J) = AIJ A(J,I) = AIJ 10 CONTINUE B(I) = T 20 CONTINUE 30 CONTINUE T = 1.E0/FLOAT(2*N-1) A(N,N) = T B(N) = T C RETURN END C COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730 COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730 C C PROGRAM: MAIN PROGRAM FOR LINPACK SGECO AND SGEFA C C VERSION: 1.0 DATE: 02/22/85 C MICROSOFT FORTRAN C C DESCRIPTION: C C COMPUTES RCOND USING SGECO FOR A SET UP AS A HILBERT C MATRIX OF SPECIFIED ORDER USING SUBROUTINE SYSGEN. C C AUTHOR: C C ADAM FRITZ C 133 MAIN STREET C AFTON, NEW YORK 13730 C C REFERENCE: C C LINPACK USERS' GUIDE C J.J. DONGARRA, C.B. MOLER, J.R. BUNCH, AND G.W. STEWART C SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS C 1979 C C SUBPROGRAMS C C SYSGEN C SGECO (FROM LINPACK) C SGEFA (FROM LINPACK) C OUT (FROM SICE) C ISAMAX (FROM BLAS) C SASUM (FROM BLAS) C SAXPY (FROM BLAS) C SDOT (FROM BLAS) C SSCAL (FROM BLAS) C SSWAP (FROM BLAS) C C VARIABLES: C C RCOND - CONDITION MEASURE C C LDA - LEADING DIMENSION OF A C N - ORDER OF A C A - INPUT: REAL SINGLE PRECISION MATRIX C OUTPUT: PACKED LU FACTORIZATION (SEE REFERENCE) C B - INPUT: REAL SINGLE PRECISION VECTOR C OUTPUT: SOLUTION VECTOR (X) C C IPVT - PIVOT INDEX VECTOR C WORK - WORK VECTOR C C INFO - STATUS VARIABLE FROM SGEFA C JOB - CONTROL VARIABLE TO SGESL C C LIN - FORTRAN LOGICAL CONSOLE INPUT C LOUT - FORTRAN LOGICAL CONSOLE OUTPUT C PRINT - PRINT CONTROL FLAG C INTEGER LDA, N, IPVT(10), INFO, JOB, LOUT, PRINT REAL A(10,10), B(10), WORK(10), RCOND, RCOND1 C LDA = 10 LIN = 3 LOUT = 3 ISUM = 0 C WRITE (LOUT, 111) 111 FORMAT (45H ***** TEST PROGRAM FOR SGECO AND SGEFA *****/) C C GET ORDER AND PRINT CODE C 10 CONTINUE WRITE (LOUT, 113) 113 FORMAT (8H ORDER: ) READ (LIN, 115) N 115 FORMAT (I2) IF (N .LT. 1 .OR. N .GT. LDA) GO TO 10 C 20 CONTINUE WRITE (LOUT, 117) 117 FORMAT (13H PRINT CODE: ) READ (LIN, 119) PRINT 119 FORMAT (I1) IF (PRINT .LT. 0 .OR. PRINT .GT. 1) GO TO 20 C C GENERATE TEST SYSTEM C CALL SYSGEN(A, LDA, N, B) IF (PRINT .EQ. 0) GO TO 30 WRITE (LOUT, 121) 121 FORMAT (28H ORIGINAL SYSTEM (BY COLUMN)/) CALL OUT(A, LDA, N, N) CALL OUT(B, LDA, N, 1) 30 CONTINUE C C COMPUTE RCOND AND FOLDED RCOND C CALL SGECO(A, LDA, N, IPVT, RCOND, WORK) IF (PRINT .EQ. 0) GO TO 40 RCOND1 = (1.0E0 + RCOND) - 1.0E0 WRITE (LOUT, 123) RCOND, RCOND1 123 FORMAT (8H RCOND: ,E16.8,9H, RCOND: ,E16.8/) 40 CONTINUE C WRITE (LOUT, 125) 125 FORMAT (24H ***** END OF TEST *****//) C STOP END C COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730  SUBROUTINE SGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(1),JOB REAL A(LDA,1),B(1) C C SGESL SOLVES THE REAL SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY SGECO OR SGEFA. C C ON ENTRY C C A REAL(LDA, N) C THE OUTPUT FROM SGECO OR SGEFA, C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER(N) C THE PIVOT VECTOR FROM SGECO OR SGEFA. C C B REAL(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF SGECO HAS SET RCOND .GT. 0.0 C OR SGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL SGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL SGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SDOT C C INTERNAL VARIABLES C REAL SDOT,T INTEGER K,KB,L,NM1 C NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL SAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL SAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = SDOT(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + SDOT(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730 C C PROGRAM: MAIN PROGRAM FOR LINPACK SGEFA AND SGESL TEST C C VERSION: 1.0 DATE: 02/22/85 C MICROSOFT FORTRAN C C DESCRIPTION: C C SOLVES; C C A * X = B C C WHERE; C C A IS SET UP AS A HILBERT MATRIX AND B IS A 'HILBERT' C VECTOR OF SPECIFIED ORDER USING SUBROUTINE SYSGEN, C X IS COMPUTED USING SGEFA TO FACTOR A, SGESL SOLVES C THE LINEAR SYSTEM IN SINGLE PRECISION, AND RESULTS C ARE SHOWN CONDITIONED ON A PRINT FLAG. C C AUTHOR: C C ADAM FRITZ C 133 MAIN STREET C AFTON, NEW YORK 13730 C C REFERENCE: C C LINPACK USERS' GUIDE C J.J. DONGARRA, C.B. MOLER, J.R. BUNCH, AND G.W. STEWART C SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS C 1979 C C SUBPROGRAMS C C SYSGEN C SGEFA (FROM LINPACK) C SGESL (FROM LINPACK) C OUT (FROM SICE) C ISAMAX (FROM BLAS) C SAXPY (FROM BLAS) C SDOT (FROM BLAS) C SSCAL (FROM BLAS) C C VARIABLES: C C LDA - LEADING DIMENSION OF A C SET TO 75 TO ASSURE LOAD IN 64K BYTES. C N - ORDER OF A C A - INPUT: REAL SINGLE PRECISION MATRIX C OUTPUT: PACKED LU FACTORIZATION (SEE REFERENCE) C B - INPUT: REAL SINGLE PRECISION VECTOR C OUTPUT: SOLUTION VECTOR (X) C C IPVT - PIVOT INDEX VECTOR C C INFO - STATUS VARIABLE FROM SGEFA C JOB - CONTROL VARIABLE TO SGESL C C LIN - FORTRAN LOGICAL CONSOLE INPUT C LOUT - FORTRAN LOGICAL CONSOLE OUTPUT C PRINT - PRINT CONTROL FLAG C INTEGER LDA, N, IPVT(75), INFO, JOB, LOUT, PRINT REAL A(75,75), B(75) C LDA = 75 LIN = 3 LOUT = 3 ISUM = 0 C WRITE (LOUT, 111) 111 FORMAT (45H ***** TEST PROGRAM FOR SGEFA AND SGESL *****/) C C GET ORDER AND PRINT CODE C 10 CONTINUE WRITE (LOUT, 113) 113 FORMAT (8H ORDER: ) READ (LIN, 115) N 115 FORMAT (I2) IF (N .LT. 1 .OR. N .GT. LDA) GO TO 10 C 20 CONTINUE WRITE (LOUT, 117) 117 FORMAT (13H PRINT CODE: ) READ (LIN, 119) PRINT 119 FORMAT (I1) IF (PRINT .LT. 0 .OR. PRINT .GT. 1) GO TO 20 C C GENERATE TEST SYSTEM C CALL SYSGEN(A, LDA, N, B) IF (PRINT .EQ. 0) GO TO 30 WRITE (LOUT, 121) 121 FORMAT (28H ORIGINAL SYSTEM (BY COLUMN)/) CALL OUT(A, LDA, N, N) CALL OUT(B, LDA, N, 1) 30 CONTINUE C C FACTOR THE MATRIX C CALL SGEFA(A, LDA, N, IPVT, INFO) IF (PRINT .EQ. 0) GO TO 40 WRITE (LOUT, 123) 123 FORMAT (28H FACTORED MATRIX (BY COLUMN)/) CALL OUT(A, LDA, N, N) 40 CONTINUE C C SOLVE THE LINEAR SYSTEM C IF (INFO .NE. 0) GO TO 60 CALL SGESL(A, LDA, N, IPVT, B, 0) IF (PRINT .EQ. 0) GO TO 50 WRITE (LOUT,125) 125 FORMAT (16H SOLUTION VECTOR/) CALL OUT(B, LDA, N, 1) 50 CONTINUE GO TO 70 60 CONTINUE C C ZERO DIAGONAL ELEMENT PRODUCED C WRITE (LOUT, 127) 127 FORMAT (25H POSSIBLE SINGULAR SYSTEM/) C 70 CONTINUE WRITE (LOUT, 129) 129 FORMAT (24H ***** END OF TEST *****//) C STOP END C COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730 COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730 C C PROGRAM: MAIN PROGRAM FOR LINPACK SGECO AND SGEDI TEST C C VERSION: 1.0 DATE: 02/22/85 C MICROSOFT FORTRAN C C DESCRIPTION: C C USES LINPACK SGECO TO COMPUTE RCOND FOR MATRIX C A AND LINPACK SGEDI TO COMPUTE DETERMINANT AND C INVERSE OF MATRIX A WHERE A IS SET UP AS A C HILBERT MATRIX OF SPECIFIED ORDER USING SYSGEN. C USES SGEFA AND SGEDI TO INVERT INVERSE TO COMPARE C WITH ORIGINAL. C C AUTHOR: C C ADAM FRITZ C 133 MAIN STREET C AFTON, NEW YORK 13730 C C REFERENCE C C LINPACK USERS' GUIDE C J.J. DONGARRA, C.B. MOLER, J.R. BUNCH, AND G.W. STEWART C SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS C 1979 C C SUBPROGRAMS C C SYSGEN C SGECO (FROM LINPACK) C SGEFA (FROM LINPACK) C SGEDI (FROM LINPACK) C OUT (FROM SICE) C ISAMAX (FROM BLAS) C SAXPY (FROM BLAS) C SDOT (FROM BLAS) C SCOPY (FROM BLAS) C SSCAL (FROM BLAS) C C VARIABLES: C C LDA - LEADING DIMENSION OF A C N - ORDER OF A C A - INPUT: REAL SINGLE PRECISION MATRIX C OUTPUT: PACKED LU FACTORIZATION (SEE REFERENCE) C A1 - COPY OF ORIGINAL MATRIX C P - PRODUCT OF ORIGINAL AND INVERSE C B - INPUT: REAL SINGLE PRECISION VECTOR C OUTPUT: SOLUTION VECTOR (X) C WORK - SCRATCH VECTOR C RCOND - MATRIX A CONDITION ESTIMATE C RCOND1 - FOLDED CONDITION ESTIMATE C C DET - DETERMINANT = DET(1) * 10**DET(2) C C IPVT - PIVOT INDEX VECTOR C C INFO - STATUS VARIABLE FROM SGEFA C JOB - CONTROL VARIABLE TO SGESL C C LIN - FORTRAN LOGICAL CONSOLE INPUT C LOUT - FORTRAN LOGICAL CONSOLE OUTPUT C PRINT - PRINT CONTROL FLAG C INTEGER I, J, LDA, N, IPVT(10), INFO, JOB, LOUT, PRINT REAL RCOND, RCOND1, DET(2) REAL A(10,10), A1(10,10), B(10), P(10,10), WORK(10) C LDA = 10 LIN = 3 LOUT = 3 ISUM = 0 C WRITE (LOUT, 111) 111 FORMAT (45H ***** TEST PROGRAM FOR SGECO AND SGEDI *****/) C C GET ORDER AND PRINT CODE C 10 CONTINUE WRITE (LOUT, 113) 113 FORMAT (8H ORDER: ) READ (LIN, 115) N 115 FORMAT (I2) IF (N .LT. 1 .OR. N .GT. LDA) GO TO 10 C 20 CONTINUE WRITE (LOUT, 117) 117 FORMAT (13H PRINT CODE: ) READ (LIN, 119) PRINT 119 FORMAT (I1) IF (PRINT .LT. 0 .OR. PRINT .GT. 1) GO TO 20 C C GENERATE TEST SYSTEM C CALL SYSGEN(A, LDA, N, B) IF (PRINT .EQ. 0) GO TO 30 WRITE (LOUT, 121) 121 FORMAT (28H ORIGINAL SYSTEM (BY COLUMN)/) CALL OUT(A, LDA, N, N) CALL OUT(B, LDA, N, 1) 30 CONTINUE C C COPY ORIGINAL MATRIX BY COLUMNS C DO 40 I = 1, N CALL SCOPY(N, A(1,I), 1, A1(1,I), 1) 40 CONTINUE C C COMPUTE CONDITION C CALL SGECO(A, LDA, N, IPVT, RCOND, WORK) IF (PRINT .EQ. 0) GO TO 50 RCOND1 = (1.0E0 + RCOND) - 1.0E0 WRITE (LOUT, 123) RCOND, RCOND1 123 FORMAT (8H RCOND: ,E16.8,10H , RCOND: ,E16.8/) 50 CONTINUE C C COMPUTE DETERMINANT AND INVERSE C IF (RCOND .EQ. 0.0E0) GO TO 120 JOB = 11 CALL SGEDI(A, LDA, N, IPVT, DET, WORK, JOB) IF (PRINT .EQ. 0) GO TO 60 WRITE (LOUT,125) 125 FORMAT (20H INVERSE (BY COLUMN)/) CALL OUT(A, LDA, N, N) WRITE (LOUT, 127) DET(1), DET(2) 127 FORMAT (14H DETERMINANT: ,F12.8,1HE,F3.0/) 60 CONTINUE C C ORIGINAL TIMES INVERSE C DO 70 I = 1, N DO 70 J = 1, N P(I,J) = SDOT(N, A1(I,1), LDA, A(1,J), 1) 70 CONTINUE IF (PRINT .EQ. 0) GO TO 80 WRITE (LOUT, 129) 129 FORMAT (35H ORIGINAL TIMES INVERSE (BY COLUMN)/) CALL OUT(P, LDA, N, N) 80 CONTINUE C C TRY TO RESTORE ORIGINAL C CALL SGEFA(A, LDA, N, IPVT, INFO) IF (INFO .GT. 0) GO TO 100 CALL SGEDI(A, LDA, N, IPVT, DET, WORK, JOB) IF (PRINT .EQ. 0) GO TO 90 WRITE (LOUT, 131) 131 FORMAT (31H INVERSE OF INVERSE (BY COLUMN)/) CALL OUT(A, LDA, N, N) 90 CONTINUE GO TO 110 100 CONTINUE WRITE (LOUT, 133) 133 FORMAT (28H ERROR: INVERSE IS SINGULAR./) 110 CONTINUE GO TO 130 C C ZERO DIAGONAL ELEMENT PRODUCED C 120 CONTINUE WRITE (LOUT, 135) 135 FORMAT (25H POSSIBLE SINGULAR SYSTEM/) C 130 CONTINUE WRITE (LOUT, 137) 137 FORMAT (24H ***** END OF TEST *****//) C STOP END C COPYRIGHT (C) 1985 ADAM FRITZ, 133 MAIN ST., AFTON, NY 13730  SUBROUTINE SGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(1),JOB REAL A(LDA,1),DET(2),WORK(1) C C SGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX C USING THE FACTORS COMPUTED BY SGECO OR SGEFA. C C ON ENTRY C C A REAL(LDA, N) C THE OUTPUT FROM SGECO OR SGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM SGECO OR SGEFA. C C WORK REAL(N) C WORK VECTOR. CONTENTS DESTROYED. C C JOB INTEGER C = 11 BOTH DETERMINANT AND INVERSE. C = 01 INVERSE ONLY. C = 10 DETERMINANT ONLY. C C ON RETURN C C A INVERSE OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE UNCHANGED. C C DET REAL(2) C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED. C OTHERWISE NOT REFERENCED. C DETERMINANT = DET(1) * 10.0**DET(2) C WITH 1.0 .LE. ABS(DET(1)) .LT. 10.0 C OR DET(1) .EQ. 0.0 . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED. C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY C AND OF SGECO HAS SET RCOND .GT. 0.0 OR SGEFA HAS SET C INFO .EQ. 0 . C C LINPACK. THIS VERSION DATED 08/14/78 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SSCAL,SSWAP C FORTRAN ABS,MOD C C INTERNAL VARIABLES C REAL T REAL TEN INTEGER I,J,K,KB,KP1,L,NM1 C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0E0 DET(2) = 0.0E0 TEN = 10.E0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = A(I,I)*DET(1) C ...EXIT IF (DET(1) .EQ. 0.0E0) GO TO 60 10 IF (ABS(DET(1)) .GE. 1.0E0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0E0 GO TO 10 20 CONTINUE 30 IF (ABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0E0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE C C COMPUTE INVERSE(U) C IF (MOD(JOB,10) .EQ. 0) GO TO 150 DO 100 K = 1, N A(K,K) = 1.0E0/A(K,K) T = -A(K,K) CALL SSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0E0 CALL SAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE C C FORM INVERSE(U)*INVERSE(L) C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 140 DO 130 KB = 1, NM1 K = N - KB KP1 = K + 1 DO 110 I = KP1, N WORK(I) = A(I,K) A(I,K) = 0.0E0 110 CONTINUE DO 120 J = KP1, N T = WORK(J) CALL SAXPY(N,T,A(1,J),1,A(1,K),1) 120 CONTINUE L = IPVT(K) IF (L .NE. K) CALL SSWAP(N,A(1,K),1,A(1,L),1) 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END C BLAS - BASIC LINEAR ALGEBRA SUBPROGRAMS C C ISAMAX C SASUM C SAXPY C SCOPY C SDOT C SNRM2 C SROT C SROTG C SSCAL C SSWAP C C REFERENCE C C LINPACK USERS' GUIDE C J.J. DONGARRA, C.B. MOLER, J.R. BUNCH, G.W. STEWART C SOCIETY FOR INDUSTRIAL AND APPLIED MATHEMATICS C 1979 C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C INTEGER FUNCTION ISAMAX(N, SX, INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SMAX INTEGER I, INCX, IX, N C ISAMAX = 0 IF (N.LT.1) RETURN ISAMAX = 1 IF (N.EQ.1) RETURN IF (INCX.EQ.1) GO TO 30 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = ABS(SX(1)) IX = IX + INCX DO 20 I=2,N IF (ABS(SX(IX)).LE.SMAX) GO TO 10 ISAMAX = I SMAX = ABS(SX(IX)) 10 IX = IX + INCX 20 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 30 SMAX = ABS(SX(1)) DO 40 I=2,N IF (ABS(SX(I)).LE.SMAX) GO TO 40 ISAMAX = I SMAX = ABS(SX(I)) 40 CONTINUE RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C REAL FUNCTION SASUM(N, SX, INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), STEMP INTEGER I, INCX, M, MP1, N, NINCX C SASUM = 0.0E0 STEMP = 0.0E0 IF (N.LE.0) RETURN IF (INCX.EQ.1) GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I=1,NINCX,INCX STEMP = STEMP + ABS(SX(I)) 10 CONTINUE SASUM = STEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,6) IF (M.EQ.0) GO TO 40 DO 30 I=1,M STEMP = STEMP + ABS(SX(I)) 30 CONTINUE IF (N.LT.6) GO TO 60 40 MP1 = M + 1 DO 50 I=MP1,N,6 STEMP = STEMP + ABS(SX(I)) + ABS(SX(I+1)) + ABS(SX(I+2)) + * ABS(SX(I+3)) + ABS(SX(I+4)) + ABS(SX(I+5)) 50 CONTINUE 60 SASUM = STEMP RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOP FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(1), SA INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N.LE.0) RETURN IF (SA.EQ.0.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I=1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF (M.EQ.0) GO TO 40 DO 30 I=1,M SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF (N.LT.4) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I+1) = SY(I+1) + SA*SX(I+1) SY(I+2) = SY(I+2) + SA*SX(I+2) SY(I+3) = SY(I+3) + SA*SX(I+3) 50 CONTINUE RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C SUBROUTINE SCOPY(N, SX, INCX, SY, INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO 1. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(1) INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I=1,N SY(IY) = SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF (M.EQ.0) GO TO 40 DO 30 I=1,M SY(I) = SX(I) 30 CONTINUE IF (N.LT.7) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,7 SY(I) = SX(I) SY(I+1) = SX(I+1) SY(I+2) = SX(I+2) SY(I+3) = SX(I+3) SY(I+4) = SX(I+4) SY(I+5) = SX(I+5) SY(I+6) = SX(I+6) 50 CONTINUE RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C REAL FUNCTION SDOT(N, SX, INCX, SY, INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(1), STEMP INTEGER I, INCX, INCY, IX, IY, M, MP1, N C STEMP = 0.0E0 SDOT = 0.0E0 IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I=1,N STEMP = STEMP + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE SDOT = STEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF (M.EQ.0) GO TO 40 DO 30 I=1,M STEMP = STEMP + SX(I)*SY(I) 30 CONTINUE IF (N.LT.5) GO TO 60 40 MP1 = M + 1 DO 50 I=MP1,N,5 STEMP = STEMP + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2) * + SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4) 50 CONTINUE 60 SDOT = STEMP RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C SUBROUTINE SROT(N, SX, INCX, SY, INCY, C, S) C C APPLIES A PLANE ROTATION. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(1), STEMP, C, S INTEGER I, INCX, INCY, IX, IY, N C IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I=1,N STEMP = C*SX(IX) + S*SY(IY) SY(IY) = C*SY(IY) - S*SX(IX) SX(IX) = STEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I=1,N STEMP = C*SX(I) + S*SY(I) SY(I) = C*SY(I) - S*SX(I) SX(I) = STEMP 30 CONTINUE RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C SUBROUTINE SROTG(SA, SB, C, S) C C CONSTRUCT GIVENS PLANE ROTATION. C JACK DONGARRA, LINPACK, 3/11/78, MODIFIED 6/82. C REAL SA, SB, C, S, ROE, SCALE, R, Z C ROE = SB IF (ABS(SA).GT.ABS(SB)) ROE = SA SCALE = ABS(SA) + ABS(SB) IF (SCALE.NE.0.0) GO TO 10 C = 1.0 S = 0.0 R = 0.0 GO TO 20 10 R = SCALE*SQRT((SA/SCALE)**2+(SB/SCALE)**2) R = SIGN(1.0,ROE)*R C = SA/R S = SB/R 20 Z = 1.0 IF (ABS(SA).GT.ABS(SB) .OR. SA*SB.EQ.0.0) Z = S IF (ABS(SB).GE.ABS(SA) .AND. ABS(SA).GT.0.0) Z = 1.0/C SA = R SB = Z RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C SUBROUTINE SSCAL(N, SA, SX, INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SA, SX(1) INTEGER I, INCX, M, MP1, N, NINCX C IF (N.LE.0) RETURN IF (INCX.EQ.1) GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I=1,NINCX,INCX SX(I) = SA*SX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF (M.EQ.0) GO TO 40 DO 30 I=1,M SX(I) = SA*SX(I) 30 CONTINUE IF (N.LT.5) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,5 SX(I) = SA*SX(I) SX(I+1) = SA*SX(I+1) SX(I+2) = SA*SX(I+2) SX(I+3) = SA*SX(I+3) SX(I+4) = SA*SX(I+4) 50 CONTINUE RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C SUBROUTINE SSWAP(N, SX, INCX, SY, INCY) C C INTERCHANGES TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO 1. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SX(1), SY(1), STEMP INTEGER I, INCX, INCY, IX, IY, M, MP1, N C IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 10 I=1,N STEMP = SX(IX) SX(IX) = SY(IY) SY(IY) = STEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,3) IF (M.EQ.0) GO TO 40 DO 30 I=1,M STEMP = SX(I) SX(I) = SY(I) SY(I) = STEMP 30 CONTINUE IF (N.LT.3) RETURN 40 MP1 = M + 1 DO 50 I=MP1,N,3 STEMP = SX(I) SX(I) = SY(I) SY(I) = STEMP STEMP = SX(I+1) SX(I+1) = SY(I+1) SY(I+1) = STEMP STEMP = SX(I+2) SX(I+2) = SY(I+2) SY(I+2) = STEMP 50 CONTINUE RETURN END C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  SUBROUTINE OUT(A, LDA, N, M) INTEGER I, J, K, M, N, IC, LDA, ICB, ICE REAL A(LDA,1) C C GENERAL PURPOSE OUTPUT ROUTINE FOR THE MATRIX A WHICH C HAS LEADING DIMENSION LDA AND IS N BY M. C IF (M.EQ.1) GO TO 30 IC = (M+4)/5 ICB = 1 ICE = 5 DO 20 K=1,IC IF (ICE.GT.M) ICE = M DO 10 I=1,N WRITE (*,99999) (A(I,J),J=ICB,ICE) 99999 FORMAT (1P5E16.8) 10 CONTINUE ICB = ICB + 5 ICE = ICE + 5 WRITE (*,99998) 99998 FORMAT (/) 20 CONTINUE WRITE (*,99997) 99997 FORMAT (/) RETURN 30 CONTINUE WRITE (*,99999) (A(I,1),I=1,N) WRITE (*,99997) RETURN END