C ALGORITHM 674, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 15, NO. 2, PP. 381-396. README.DOC This disk accompanies the submission "Fortran codes for estimating the one- norm of a real or complex matrix, with applications to condition estimation" to ACM Trans. Math. Soft. The files are: SONRUN.FOR Test driver program for subroutine SONEST. RAND.FOR Random number generators. SONEST.FOR Subroutine SONEST. LPSTUFF.FOR LINPACK SGE?? routines (for real matrices) and associated BLAS routines. CONRUN.FOR Test driver program for subroutine CONEST. CONEST.FOR Subroutine CONEST. LPCMPLX.FOR LINPACK CGE?? routines (for complex matrices) and associated BLAS routines. MYCBLAS.FOR Modified complex BLAS routines. The dependency among these files is as follows: SONRUN.FOR uses the subroutines and functions contained in RAND.FOR, SONEST.FOR and LPSTUFF.FOR CONRUN.FOR uses the subroutines and functions contained in RAND.FOR, CONEST.FOR, LPCMPLX.FOR and MYCBLAS.FOR Nicholas J. Higham na.higham@score.stanford.edu University of Manchester Manchester M13 9PL England 17th March 1988. PROGRAM CONRUN C C THIS IS A TEST DRIVER FOR CONEST. C IT APPLIES CONEST TO SEVERAL CLASSES OF MATRIX AND EVALUATES C VARIOUS PERFORMANCE MEASURES. IT SERVES THE FOLLOWING PURPOSES. C C (1) ILLUSTRATES HOW TO USE CONEST (WITH LINPACK'S CGEFA/CGESL). C C (2) RIGOROUSLY TESTS CONEST. EXERCISES ALL STATEMENTS AND C FEATURES OF THE CODE. C C (3) IS A TOOL FOR STUDYING THE BEHAVIOUR OF CONEST. C C RESULTS ARE REPORTED SEPARATELY FOR EACH OF FIVE CLASSES OF TEST C MATRIX. WITHIN EACH CLASS THE MINIMUM, MAXIMUM, AND AVERAGE VALUES C OF SEVEN STATISTICS ARE REPORTED FOR EACH COMBINATION OF CONDITION C NUMBER (WHERE APPLICABLE) AND N. C C TO CHECK THAT CONEST IS WORKING CORRECTLY CHECK THE FOLLOWING. C - STATISTIC 1 (CONEST ESTIMATE/TRUE NORM). C THE MINIMUM VALUES SHOULD GENERALLY BE AT LEAST 0.1. C THE POSSIBILITY OF SMALLER VALUES CANNOT BE RULED C OUT, SINCE THE ALGORITHM UNDERLYING THE CODE IS C NOT COMPLETELY RELIABLE. C - STATISTIC 4 (NULL VECTOR RATIO). C THE MINIMUM VALUE SHOULD BE 1.0, EXCEPT THAT THERE C WILL BE SOME DETERIORATION FOR LARGER CONDITION C NUMBERS, AS NOTED BELOW. C - THE SUMMARY OF THE NUMBER OF ITERATIONS. THIS SHOULD SHOW C RELATIVELY FEW OCCURRENCES OF 4 OR 5 ITERATIONS. C C THE OTHER STATISTICS ARE OF INTEREST FOR INVESTIGATING AND C UNDERSTANDING THE PERFORMANCE OF CONEST. FOR EXAMPLE, ONE C MIGHT EXAMINE THE FOLLOWING ASPECTS. C C - THE NUMBER OF TIMES THE EXTRA ESTIMATE IS USED. C - THE FREQUENCY OF EXACT ESTIMATES (VALUE OF 1.0 FOR STATISTIC 1). C - HOW THE ABOVE FACTORS VARY WITH N, THE CONDITION NUMBER, AND THE C TYPE OF MATRIX. C C THE USER SHOULD ADJUST THE PARAMETERS NCOND, NOFN, ITIMES, AND THE C ARRAY NVALS, TO TASTE - E.G. TO SUIT THE COMPUTER TIME AVAILABLE. C THE MODEST DEFAULT VALUES BELOW ARE SUITABLE FOR AN INITIAL TEST C RUN. FOR EFFICIENCY, THE CASES N=1, AND MODE=-2, SHOULD BE RUN C WITH ITIMES=1, SEPARATELY FROM THE MAIN TESTING WITH ITIMES>1 C C THIS VERSION DATED MARCH 16TH, 1988. C NICK HIGHAM, UNIVERSITY OF MANCHESTER. C C MACHINE DEPENDENCE C C THE ONLY MACHINE DEPENDENT PARAMETER IS EPS IN THE FOLLOWING C PARAMETER STATEMENT. SET THIS TO THE MACHINE PRECISION C (UNIT ROUNDOFF). PARAMETER (EPS=2.0**(-23)) C C SUBROUTINES AND FUNCTIONS. C C LINPACK - CGECO, CGEFA, CGESL, CGEDI. C BLAS - THOSE REQUIRED BY THE ABOVE LINPACK ROUTINES AND BY CONEST. C RANDOM NUMBER GENERATORS URAND AND RANDN, SUPPLIED WITH THIS CODE. C C NOTE THAT CGECO USES A MODIFIED VERSION OF THE 1-NORM C SO THE RATIOS REPORTED IN STATISTICS 5 AND 6 ARE NOT THOSE OF C STRICTLY COMPARABLE QUANTITIES. C C THE TEST MATRICES FOR MODE.GE.0 WILL HAVE 2-NORM CONDITION NUMBERS C RANGING FROM 1 TO CONDMX, TAKING ON NCOND DIFFERENT VALUES THAT C ARE POWERS OF CONDPW (A REAL VARIABLE SET BELOW). C THE MAXIMUM CONDITION NUMBER, CONDMX, IS JUST BELOW THE LEVEL C CORRESPONDING TO THE MATRIX BEING 'SINGULAR TO WORKING PRECISION'. C C N.B. FOR CONDITION NUMBERS ON THE ORDER OF 1.0/EPS, ROUNDING ERRORS C ARE SIGNIFICANT AND THE RATIOS NEED TO BE INTERPRETED WITH CARE. C STATISTICS 4, 5 AND 6 (NULL VECTOR RATIO, CGECO/TRUE NORM C AND CGECO/CONEST) ARE ALL 'UNRELIABLE' FOR LARGE CONDITION C NUMBERS. HOWEVER, STATISTIC 1 (CONEST/TRUE NORM) DOES NOT C DETERIORATE FOR LARGE CONDITION NUMBERS BECAUSE OF CORRELATION C WITH THE 'TRUE' CONDITION NUMBER OBTAINED VIA CGEDI (ESSENTIALLY C THESE TWO ESTIMATES USE THE SAME COMPUTED COLUMNS OF INV(A)). PARAMETER (CONDMX = (1.0/EPS)/10.0) PARAMETER (NCOND=3) C C NUMBER OF N VALUES PARAMETER (NOFN=3) C NUMBER OF RANDOM REPETITIONS FOR EACH COMBINATION OF MATRIX TYPE, C N, AND CONDITION NUMBER. PARAMETER (ITIMES=2) C C NUMBER OF STATISTICS TO INVESTIGATE PARAMETER (NESTS=7) C CONEST TAKES BETWEEN 1 AND ITMAX ITERATIONS, EXCEPT C IN THE TRIVIAL CASE N=1 WHEN IT TAKES 'HALF AN ITERATION' C (WHICH WE COUNT HERE AS ZERO ITERATIONS). PARAMETER (ITMAX=5) C C LDA MUST BE AT LEAST AS BIG AS THE MAXIMUM VALUE OF N. PARAMETER (LDA=50, LDASQ=LDA*LDA) C COMPLEX A(LDA,LDA), X(LDA), V(LDA), ASAVE(LDA,LDA) COMPLEX W1(LDA), W2(LDA), DET(2), CSUM REAL S(NESTS,ITIMES), RES(NESTS,NOFN,NCOND,3) INTEGER NVALS(NOFN), IPIVOT(LDA), ITS(0:ITMAX) C C MULTIPLICATIVE FACTOR FOR GENERATING CONDITION NUMBERS. CONDPW = CONDMX**(1.0/REAL(NCOND-1)) C C THE NOFN DIFFERENT VALUES OF N TO BE USED. DATA NVALS/1,5,10/ * C INITIALISE A TO ZERO. DATA A/LDASQ*(0.0,0.0)/ C C FOR EACH OF THE FIVE CLASSES OF MATRIX... DO 60,MODE = -2,2 C ISING = 0 IF (MODE.LT.0) THEN C FOR MODE.LT.0 THE CONDITION NUMBER IS NOT A PARAMETER. INCOND = 1 ELSE INCOND = NCOND END IF C C RESET CONEST ITERATION STATISTICS. DO 1,I=0,ITMAX ITS(I) = 0 1 CONTINUE ICYCLE = 0 IXTRA = 0 C C FOR EACH VALUE OF N... DO 30,I=1,NOFN N=NVALS(I) C FOR EACH CONDITION NUMBER... DO 30,J=1,INCOND C DO 10,K=1,ITIMES C SEED FOR THE RANDOM NUMBER GENERATOR. IRAND = K*J*I+1 3 IF (MODE.GE.0) THEN COND = CONDPW**(J-1) CALL UDV(N,A,LDA,COND,MODE,W1,W2,IRAND) ELSEIF (MODE.EQ.-1) THEN CALL RANDA (N,A,LDA,IRAND) ELSEIF (MODE.EQ.-2) THEN CALL BIDIAG (N,A,LDA) END IF CALL CCOPY(LDA*LDA,A,1,ASAVE,1) C ANORM = 0.0 DO 2,L=1,N ANORM = MAX(ANORM, SCSUM1(N,A(1,L),1)) 2 CONTINUE C GET PA=LU FACTORISATION TOGETHER WITH LINPACK CONDITION ESTIMATE. CALL CGECO(A,LDA,N,IPIVOT,RCOND,W1) C C IN THE UNLIKELY EVENT OF A SINGULAR U IN PA=LU (FOR MODE.GE.-1), C GENERATE ANOTHER A. IF (RCOND .EQ. 0) THEN ISING = ISING+1 GOTO 3 END IF C C ------------------------ C THIS SECTION OF CODE APPLIES CONEST. THE FIRST LINE, AND THE ONES C INVOLVING ISOLVE AND ESTOLD, COULD BE OMITTED IN PRACTICE; C THEY ARE USED HERE TO MONITOR THE PERFORMANCE OF CONEST. EST = 0.0 ISOLVE = 0 KASE = 0 5 CALL CONEST(N,V,X,EST,KASE) IF (KASE .NE. 0) THEN ISOLVE = ISOLVE+1 ESTOLD = EST CALL CGESL(A,LDA,N,IPIVOT,X,KASE-1) GOTO 5 END IF C ------------------------ C C MAKE DEDUCTIONS ABOUT HOW MANY ITERATIONS (ITER) WERE USED, AND C WHETHER CYCLING WAS DETECTED. IF (N .GT. 1) THEN IF (EST .NE. ESTOLD) IXTRA = IXTRA+1 ITSLVE = ISOLVE-1 ITER = INT( (ITSLVE+1)/2 ) IF ( MOD(ITSLVE,2) .EQ. 1 ) ICYCLE = ICYCLE+1 ITS(ITER) = ITS(ITER)+1 ELSE ITER = INT(ISOLVE/2) ITS(ITER) = ITS(ITER)+1 END IF C C CHECK APPROXIMATE NULL VECTOR. DO 6,L=1,N CSUM = (0.0,0.0) DO 8,LL=1,N CSUM = CSUM + ASAVE(L,LL)*V(LL) 8 CONTINUE X(L) = CSUM 6 CONTINUE XEST = SCSUM1(N,V,1)/SCSUM1(N,X,1) C C COMPUTE THE 'TRUE' VALUE OF NORM1(INV(A)). CALL CGEDI(A,LDA,N,IPIVOT,DET,W1,1) AINORM = 0.0 DO 7,L=1,N AINORM = MAX(AINORM,SCSUM1(N,A(1,L),1)) 7 CONTINUE COND1 = ANORM*AINORM C C FORM THE STATISTICS OF INTEREST. 9 S(1,K) = EST/AINORM S(2,K) = REAL(ITER) S(3,K) = REAL(ISOLVE) S(4,K) = XEST/EST S(5,K) = 1.0/(RCOND*COND1) S(6,K) = EST*ANORM*RCOND S(7,K) = COND1 C 10 CONTINUE C C RETAIN ONLY THE MINIMUM, MAXIMUM AND AVERAGE VALUES OF THE STATS. DO 20,L=1,NESTS SMIN=S(L,1) SMAX=S(L,1) SUM=0.0 DO 15,K=1,ITIMES SMIN=MIN(S(L,K),SMIN) SUM=SUM+S(L,K) SMAX=MAX(S(L,K),SMAX) 15 CONTINUE RES(L,I,J,1)=SMIN RES(L,I,J,2)=SMAX RES(L,I,J,3)=SUM/REAL(ITIMES) 20 CONTINUE 30 CONTINUE C C PRINT RESULTS WRITE(6,'(//)') WRITE(6,*) 'RESULTS FOR MODE = ',MODE,' ITIMES = ',ITIMES IF (MODE .EQ. -2) THEN WRITE(6,*) 'UPPER BIDIAGONAL, U(I,J) = 1' ELSE IF (MODE .EQ. -1) THEN WRITE(6,*) 'ELTS OF A FROM UNIFORM (-1,1) DISTRIBUTION' ELSE IF (MODE .EQ. 0) THEN WRITE(6,*) 'ONE LARGE SINGULAR VALUE' ELSE IF (MODE .EQ. 1) THEN WRITE(6,*) 'ONE SMALL SINGULAR VALUE' ELSE IF (MODE .EQ. 2) THEN WRITE(6,*) 'EXPONENTIALLY DISTRIBUTED SINGULAR VALUES' END IF WRITE(6,*) '##############################################' WRITE(6,*) WRITE(6,*) 'CONEST ITERATION BEHAVIOUR' WRITE(6,'(A18,6I6)') 'NO. ITERATIONS ', 0,1,2,3,4,5 WRITE(6,'(A18,6I6)') ' # TIMES', (ITS(I),I=0,ITMAX) WRITE(6,'(A33,I8)') ' NUMBER OF TIMES CYCLING DETECTED', ICYCLE WRITE(6,'(A31,I8)') ' NUMBER OF TIMES EXTRA EST USED', IXTRA WRITE(6,'(A25,I8)') ' TOTAL NUMBER OF MATRICES ', + NOFN*INCOND*ITIMES WRITE(6,'(A22,I8/)') ' NUMBER OF SINGULAR U =', ISING DO 50,L=1,NESTS WRITE(6,*) WRITE(6,*) '********************************************' WRITE(6,*) 'STATISTIC NO ',L IF (L .EQ. 1) THEN WRITE(6,*) 'RATIO - CONEST ESTIMATE / NORM1(INV(A))' ELSE IF (L .EQ. 2) THEN WRITE(6,*) 'NUMBER OF ITERATIONS' ELSE IF (L .EQ. 3) THEN WRITE(6,*) 'NUMBER OF LINEAR SYSTEMS SOLVED' ELSE IF (L .EQ. 4) THEN WRITE(6,*) 'CHECK ON APPROX NULL VECTOR - SHOULD BE 1.0' ELSE IF (L .EQ. 5) THEN WRITE(6,*) 'RATIO - LINPACK ESTIMATE / NORM1(INV(A))' WRITE(6,*) '***NB** LINPACK USES A MODIFIED 1-NORM' ELSE IF (L .EQ. 6) THEN WRITE(6,*) 'RATIO - CONEST ESTIMATE / LINPACK ESTIMATE' WRITE(6,*) '***NB** LINPACK USES A MODIFIED 1-NORM' ELSE IF (L .EQ. 7) THEN WRITE(6,*) '1-NORM CONDITION NUMBER' END IF WRITE(6,*) 'MIN, MAX THEN AVE' DO 40,K=1,3 WRITE(6,*) IF (MODE .LE. -1) THEN WRITE(6,'(A5,6I12)') ' N ', (NVALS(I),I=1,NOFN) ELSE WRITE(6,'(A12,6I12)') 'COND / N ',(NVALS(I),I=1,NOFN) END IF WRITE(6,*) '--------------------------------------' DO 35,J=1,INCOND IF (MODE .LE. -1) THEN WRITE(6,'(5X,6E12.4)') (RES(L,I,J,K),I=1,NOFN) ELSE WRITE(6,'(7E12.4)') CONDPW**(J-1),(RES(L,I,J,K),I=1,NOFN) END IF 35 CONTINUE 40 CONTINUE 50 CONTINUE 60 CONTINUE C END C SUBROUTINE UDV(N,A,LDA,COND2,MODE,X,D,IRAND) C GENERATES RANDOM A WITH PRE-ASSIGNED SINGULAR VALUE DISTRIBUTION. C A = UDV, WHERE U AND V ARE RANDOM COMPLEX ORTHOGONAL MATRICES C AND D IS DIAGONAL WITH THE SINGULAR VALUES ON ITS DIAGONAL. C U AND V ARE BOTH THE PRODUCT OF NHTS RANDOM HOUSEHOLDER MATRICES. PARAMETER (NHTS = 5) C IN THIS ROUTINE WE MAKE USE OF FORTRAN'S IMPLICIT REAL TO COMPLEX C CONVERSIONS. COMPLEX A(LDA,N), X(N), D(N), T DO 5,I=1,N DO 5,J=1,N A(I,J) = 0.0 5 CONTINUE IF (MODE.EQ.0) THEN A(1,1) = 1.0 DO 10,I=2,N A(I,I) = 1.0/COND2 10 CONTINUE ELSE IF (MODE.EQ.1) THEN DO 15,I=1,N-1 A(I,I) = 1.0 15 CONTINUE A(N,N) = 1.0/COND2 ELSE IF (MODE.EQ.2) THEN A(1,1) = 1.0 IF (N.GT.1) THEN ALPHA = COND2**(-1.0/REAL(N-1)) DO 20,I=2,N A(I,I) = ALPHA*A(I-1,I-1) 20 CONTINUE END IF END IF DO 25,I=1,NHTS CALL UTIMEA(N,A,LDA,X,D,IRAND) 25 CONTINUE C TRANSPOSE A (DON'T NEED TO CONJUGATE) C (AVOIDS NEED FOR AN ANALOGUE OF UTIMEA FOR POST-MULTIPLICATION). DO 30,I=1,N-1 DO 30,J=I+1,N T = A(I,J) A(I,J) = A(J,I) A(J,I) = T 30 CONTINUE DO 35,I=1,NHTS CALL UTIMEA(N,A,LDA,X,D,IRAND) 35 CONTINUE END C SUBROUTINE UTIMEA(N,A,LDA,V,VTA,IRAND) COMPLEX A(LDA,N), V(N), VTA(N), CSUM REAL RANDN C PRE-MULTIPLIES A BY A RANDOM COMPLEX HOUSEHOLDER MATRIX. SN2 = 0.0 DO 10,I=1,N VR = RANDN(IRAND) VI = RANDN(IRAND) V(I) = CMPLX(VR, VI) SN2 = SN2 + VR**2 + VI**2 10 CONTINUE DO 30,J=1,N CSUM = (0.0, 0.0) DO 20,I=1,N CSUM =CSUM + CONJG(V(I))*A(I,J) 20 CONTINUE VTA(J) = CSUM 30 CONTINUE DO 40,I=1,N DO 40,J=1,N A(I,J) = A(I,J) - V(I)*VTA(J) * CMPLX(2.0/SN2) 40 CONTINUE END C SUBROUTINE RANDA (N,A,LDA,IRAND) C GENERATES RANDOM COMPLEX A WITH ELEMENTS HAVING REAL AND C IMAGINARY PARTS FROM THE UNIFORM (-1,1) DISTRIBUTION. COMPLEX A(LDA,N) REAL URAND DO 10,J=1,N DO 10,I=1,N A(I,J) = CMPLX(-1.0+2.0*URAND(IRAND), -1.0+2.0*URAND(IRAND)) 10 CONTINUE END C SUBROUTINE BIDIAG (N,A,LDA) COMPLEX A(LDA,N) DO 10,J=1,N DO 10,I=1,N IF ( (J.EQ.I) .OR. (J.EQ.I+1) ) THEN A(I,J) = (1.0,0.0) ELSE A(I,J) = (0.0,0.0) END IF 10 CONTINUE END C THESE ROUTINES OBTAINED C FROM: NETLIB C DATE: FRI, 10 APR 87 08:58:02 CST C SUBJECT: CGEFA CGESL CGECO CGEDI FROM LINPACK C MODIFIED BY NICK HIGHAM TO REPLACE '1' IN LAST DIMENSIONS BY '*'. C SUBROUTINE CGECO(A,LDA,N,IPVT,RCOND,Z) INTEGER LDA,N,IPVT(*) COMPLEX A(LDA,*),Z(*) REAL RCOND C C CGECO FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, CGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B , FOLLOW CGECO BY CGESL. C TO COMPUTE INVERSE(A)*C , FOLLOW CGECO BY CGESL. C TO COMPUTE DETERMINANT(A) , FOLLOW CGECO BY CGEDI. C TO COMPUTE INVERSE(A) , FOLLOW CGECO BY CGEDI. C C ON ENTRY C C A COMPLEX(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 COMPLEX(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 CGEFA C BLAS CAXPY,CDOTC,CSSCAL,SCASUM C FORTRAN ABS,AIMAG,AMAX1,CMPLX,CONJG,REAL C C INTERNAL VARIABLES C COMPLEX CDOTC,EK,T,WK,WKM REAL ANORM,S,SCASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C COMPLEX ZDUM,ZDUM1,ZDUM2,CSIGN1 REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) CSIGN1(ZDUM1,ZDUM2) = CABS1(ZDUM1)*(ZDUM2/CABS1(ZDUM2)) C C COMPUTE 1-NORM OF A C ANORM = 0.0E0 DO 10 J = 1, N ANORM = AMAX1(ANORM,SCASUM(N,A(1,J),1)) 10 CONTINUE C C FACTOR C CALL CGEFA(A,LDA,N,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND CTRANS(A)*Y = E . C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A . C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W = E . C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE CTRANS(U)*W = E C EK = (1.0E0,0.0E0) DO 20 J = 1, N Z(J) = (0.0E0,0.0E0) 20 CONTINUE DO 100 K = 1, N IF (CABS1(Z(K)) .NE. 0.0E0) EK = CSIGN1(EK,-Z(K)) IF (CABS1(EK-Z(K)) .LE. CABS1(A(K,K))) GO TO 30 S = CABS1(A(K,K))/CABS1(EK-Z(K)) CALL CSSCAL(N,S,Z,1) EK = CMPLX(S,0.0E0)*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = CABS1(WK) SM = CABS1(WKM) IF (CABS1(A(K,K)) .EQ. 0.0E0) GO TO 40 WK = WK/CONJG(A(K,K)) WKM = WKM/CONJG(A(K,K)) GO TO 50 40 CONTINUE WK = (1.0E0,0.0E0) WKM = (1.0E0,0.0E0) 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + CABS1(Z(J)+WKM*CONJG(A(K,J))) Z(J) = Z(J) + WK*CONJG(A(K,J)) S = S + CABS1(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*CONJG(A(K,J)) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0E0/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) C C SOLVE CTRANS(L)*Y = W C DO 120 KB = 1, N K = N + 1 - KB IF (K .LT. N) Z(K) = Z(K) + CDOTC(N-K,A(K+1,K),1,Z(K+1),1) IF (CABS1(Z(K)) .LE. 1.0E0) GO TO 110 S = 1.0E0/CABS1(Z(K)) CALL CSSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0E0/SCASUM(N,Z,1) CALL CSSCAL(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 CAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (CABS1(Z(K)) .LE. 1.0E0) GO TO 130 S = 1.0E0/CABS1(Z(K)) CALL CSSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0E0/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) YNORM = S*YNORM C C SOLVE U*Z = V C DO 160 KB = 1, N K = N + 1 - KB IF (CABS1(Z(K)) .LE. CABS1(A(K,K))) GO TO 150 S = CABS1(A(K,K))/CABS1(Z(K)) CALL CSSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (CABS1(A(K,K)) .NE. 0.0E0) Z(K) = Z(K)/A(K,K) IF (CABS1(A(K,K)) .EQ. 0.0E0) Z(K) = (1.0E0,0.0E0) T = -Z(K) CALL CAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0E0/SCASUM(N,Z,1) CALL CSSCAL(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 CGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(*),JOB COMPLEX A(LDA,*),DET(2),WORK(*) C C CGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX C USING THE FACTORS COMPUTED BY CGECO OR CGEFA. C C ON ENTRY C C A COMPLEX(LDA, N) C THE OUTPUT FROM CGECO OR CGEFA. 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 CGECO OR CGEFA. C C WORK COMPLEX(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 COMPLEX(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. CABS1(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 IF CGECO HAS SET RCOND .GT. 0.0 OR CGEFA 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 CAXPY,CSCAL,CSWAP C FORTRAN ABS,AIMAG,CMPLX,MOD,REAL C C INTERNAL VARIABLES C COMPLEX T REAL TEN INTEGER I,J,K,KB,KP1,L,NM1 C COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = (1.0E0,0.0E0) DET(2) = (0.0E0,0.0E0) TEN = 10.0E0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = A(I,I)*DET(1) C ...EXIT IF (CABS1(DET(1)) .EQ. 0.0E0) GO TO 60 10 IF (CABS1(DET(1)) .GE. 1.0E0) GO TO 20 DET(1) = CMPLX(TEN,0.0E0)*DET(1) DET(2) = DET(2) - (1.0E0,0.0E0) GO TO 10 20 CONTINUE 30 IF (CABS1(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/CMPLX(TEN,0.0E0) DET(2) = DET(2) + (1.0E0,0.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,0.0E0)/A(K,K) T = -A(K,K) CALL CSCAL(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,0.0E0) CALL CAXPY(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,0.0E0) 110 CONTINUE DO 120 J = KP1, N T = WORK(J) CALL CAXPY(N,T,A(1,J),1,A(1,K),1) 120 CONTINUE L = IPVT(K) IF (L .NE. K) CALL CSWAP(N,A(1,K),1,A(1,L),1) 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END SUBROUTINE CGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(*),INFO COMPLEX A(LDA,*) C C CGEFA FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION. C C CGEFA IS USUALLY CALLED BY CGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR CGECO) = (1 + 9/N)*(TIME FOR CGEFA) . C C ON ENTRY C C A COMPLEX(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 CGESL OR CGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN CGECO 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 CAXPY,CSCAL,ICAMAX C FORTRAN ABS,AIMAG,REAL C C INTERNAL VARIABLES C COMPLEX T INTEGER ICAMAX,J,K,KP1,L,NM1 C COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) 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 = ICAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (CABS1(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,0.0E0)/A(K,K) CALL CSCAL(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 CAXPY(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 (CABS1(A(N,N)) .EQ. 0.0E0) INFO = N RETURN END SUBROUTINE CGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(*),JOB COMPLEX A(LDA,*),B(*) C C CGESL SOLVES THE COMPLEX SYSTEM C A * X = B OR CTRANS(A) * X = B C USING THE FACTORS COMPUTED BY CGECO OR CGEFA. C C ON ENTRY C C A COMPLEX(LDA, N) C THE OUTPUT FROM CGECO OR CGEFA. 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 CGECO OR CGEFA. C C B COMPLEX(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE CTRANS(A)*X = B WHERE C CTRANS(A) IS THE CONJUGATE 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 NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF CGECO HAS SET RCOND .GT. 0.0 C OR CGEFA HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL CGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1, P C CALL CGESL(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 CAXPY,CDOTC C FORTRAN CONJG C C INTERNAL VARIABLES C COMPLEX CDOTC,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 CAXPY(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 CAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE CTRANS(A) * X = B C FIRST SOLVE CTRANS(U)*Y = B C DO 60 K = 1, N T = CDOTC(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/CONJG(A(K,K)) 60 CONTINUE C C NOW SOLVE CTRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + CDOTC(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 SUBROUTINE CSCAL(N,CA,CX,INCX) C C SCALES A VECTOR BY A CONSTANT. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CA,CX(*) INTEGER I,INCX,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 CX(I) = CA*CX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N CX(I) = CA*CX(I) 30 CONTINUE RETURN END SUBROUTINE CSSCAL(N,SA,CX,INCX) C C SCALES A COMPLEX VECTOR BY A REAL CONSTANT. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*) REAL SA INTEGER I,INCX,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 CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) 30 CONTINUE RETURN END SUBROUTINE CSWAP (N,CX,INCX,CY,INCY) C C INTERCHANGES TWO VECTORS. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*),CY(*),CTEMP 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 CTEMP = CX(IX) CX(IX) = CY(IY) CY(IY) = CTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 20 DO 30 I = 1,N CTEMP = CX(I) CX(I) = CY(I) CY(I) = CTEMP 30 CONTINUE RETURN END INTEGER FUNCTION ICAMAX(N,CX,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*) REAL SMAX INTEGER I,INCX,IX,N COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C ICAMAX = 0 IF( N .LT. 1 ) RETURN ICAMAX = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = CABS1(CX(1)) IX = IX + INCX DO 10 I = 2,N IF(CABS1(CX(IX)).LE.SMAX) GO TO 5 ICAMAX = I SMAX = CABS1(CX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 SMAX = CABS1(CX(1)) DO 30 I = 2,N IF(CABS1(CX(I)).LE.SMAX) GO TO 30 ICAMAX = I SMAX = CABS1(CX(I)) 30 CONTINUE RETURN END REAL FUNCTION SCASUM(N,CX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND C RETURNS A SINGLE PRECISION RESULT. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*) REAL STEMP INTEGER I,INCX,N,NINCX C SCASUM = 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(REAL(CX(I))) + ABS(AIMAG(CX(I))) 10 CONTINUE SCASUM = STEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) 30 CONTINUE SCASUM = STEMP RETURN END SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*),CY(*),CA INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.0)RETURN IF (ABS(REAL(CA)) + ABS(AIMAG(CA)) .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 CY(IY) = CY(IY) + CA*CX(IX) 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 CY(I) = CY(I) + CA*CX(I) 30 CONTINUE RETURN END COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS, CONJUGATING THE FIRST C VECTOR. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*),CY(*),CTEMP INTEGER I,INCX,INCY,IX,IY,N C CTEMP = (0.0,0.0) CDOTC = (0.0,0.0) 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 CTEMP = CTEMP + CONJG(CX(IX))*CY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE CDOTC = CTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I = 1,N CTEMP = CTEMP + CONJG(CX(I))*CY(I) 30 CONTINUE CDOTC = CTEMP RETURN END C THESE ROUTINES OBTAINED C FROM: NETLIB C DATE: FRI, 10 APR 87 08:59:15 CST C SUBJECT: SGEFA SGECO SGESL SGEDI FROM LINPACK C MODIFIED BY NICK HIGHAM TO REPLACE '1' IN LAST DIMENSIONS BY '*'. C SUBROUTINE SGECO(A,LDA,N,IPVT,RCOND,Z) INTEGER LDA,N,IPVT(*) REAL A(LDA,*),Z(*) 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 A*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 SGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(*),JOB REAL A(LDA,*),DET(2),WORK(*) 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 IF 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 C COMPUTE DETERMINANT C IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0E0 DET(2) = 0.0E0 TEN = 10.0E0 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 SUBROUTINE SGEFA(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(*),INFO REAL A(LDA,*) 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 SUBROUTINE SGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(*),JOB REAL A(LDA,*),B(*) 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 C THE ORDER OF THE MATRIX A . C C IPVT 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 NOT 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 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(*) 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 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(*),SY(*),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 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(*),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 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = ABS(SX(1)) IX = IX + INCX DO 10 I = 2,N IF(ABS(SX(IX)).LE.SMAX) GO TO 5 ISAMAX = I SMAX = ABS(SX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 SMAX = ABS(SX(1)) DO 30 I = 2,N IF(ABS(SX(I)).LE.SMAX) GO TO 30 ISAMAX = I SMAX = ABS(SX(I)) 30 CONTINUE RETURN END 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(*),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 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(*),SY(*),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 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(*),SY(*),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 NEXT BLAS IS NOT FROM ABOVE NETLIB MESSAGE. IT'S ONE I MODIFIED C BY DELETING COMMENTS AND CHANGING TO *. * SUBROUTINE SCOPY(N,SX,INCX,SY,INCY) REAL SX(*),SY(*) INTEGER I,INCX,INCY,IX,IY,M,MP1,N IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 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 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 INTEGER FUNCTION ICMAX1(N,CX,INCX) C C FINDS THE INDEX OF ELEMENT WHOSE REAL PART HAS MAX. ABS. VALUE. C BASED ON ICAMAX BY JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED BY NICK HIGHAM, MAY 12, 1987. C COMPLEX CX(*) REAL SMAX INTEGER I,INCX,IX,N COMPLEX ZDUM REAL CABS1 C C ... NEXT LINE IS THE ONLY MODIFICATION. CABS1(ZDUM) = ABS(REAL(ZDUM)) C ICMAX1 = 0 IF( N .LT. 1 ) RETURN ICMAX1 = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = CABS1(CX(1)) IX = IX + INCX DO 10 I = 2,N IF(CABS1(CX(IX)).LE.SMAX) GO TO 5 ICMAX1 = I SMAX = CABS1(CX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 SMAX = CABS1(CX(1)) DO 30 I = 2,N IF(CABS1(CX(I)).LE.SMAX) GO TO 30 ICMAX1 = I SMAX = CABS1(CX(I)) 30 CONTINUE RETURN END C REAL FUNCTION SCSUM1(N,CX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND C RETURNS A SINGLE PRECISION RESULT. C BASED ON SCASUM BY JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED BY NICK HIGHAM, MAY 12, 1987. C THE CHANGE IS TO USE THE 'GENUINE' ABSOLUTE VALUE. C COMPLEX CX(*) REAL STEMP INTEGER I,INCX,N,NINCX C SCSUM1 = 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 C ... NEXT LINE MODIFIED. STEMP = STEMP + ABS(CX(I)) 10 CONTINUE SCSUM1 = STEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N C ... NEXT LINE MODIFIED. STEMP = STEMP + ABS(CX(I)) 30 CONTINUE SCSUM1 = STEMP RETURN END C SUBROUTINE CCOPY(N,CX,INCX,CY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C JACK DONGARRA, LINPACK, 3/11/78. C COMPLEX CX(*),CY(*) 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 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 CY(IY) = CX(IX) 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 CY(I) = CX(I) 30 CONTINUE RETURN END REAL FUNCTION URAND(IY) INTEGER IY C C URAND FROM FORSYTHE, MALCOLM AND MOLER (AS PROVIDED BY NETLIB). C C URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2. THE INTEGER IY C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL C TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED C IN THE INTERVAL (0,1). C INTEGER IA,IC,ITWO,M2,M,MIC DOUBLE PRECISION HALFM REAL S DOUBLE PRECISION DATAN,DSQRT DATA M2/0/,ITWO/2/ IF (M2 .NE. 0) GO TO 20 C C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH C M = 1 10 M2 = M M = ITWO*M2 IF (M .GT. M2) GO TO 10 HALFM = M2 C C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD C IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5 IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1 MIC = (M2 - IC) + M2 C C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT C S = 0.5/HALFM C C COMPUTE NEXT RANDOM NUMBER C 20 IY = IY*IA C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW C INTEGER OVERFLOW ON ADDITION C IF (IY .GT. MIC) IY = (IY - M2) - M2 C IY = IY + IC C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE C WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION C IF (IY/2 .GT. M2) IY = (IY - M2) - M2 C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER C OVERFLOW AFFECTS THE SIGN BIT C IF (IY .LT. 0) IY = (IY + M2) + M2 URAND = FLOAT(IY)*S RETURN END C C REAL FUNCTION RANDN(IY) C POLAR METHOD FOR NORMAL (0,1) DEVIATES. C FROM D.E. KNUTH, SEMINUMERICAL ALGORITHMS, SECOND EDITION, C ADDISON-WESLEY, 1981, PAGE 117. REAL URAND SAVE IOLD, RNDOLD DATA IOLD/0/ IF (IOLD .EQ. 0) THEN 10 V1 = 2.0*URAND(IY)-1.0 V2 = 2.0*URAND(IY)-1.0 R = V1**2 + V2**2 IF (R .GE. 1.0) GO TO 10 FACT = SQRT(-2.0*LOG(R)/R) RNDOLD = V1*FACT RANDN = V2*FACT IOLD = 1 ELSE RANDN = RNDOLD IOLD = 0 ENDIF END SUBROUTINE SONEST (N, V, X, ISGN, EST, KASE) INTEGER N, ISGN(N), KASE REAL V(N), X(N), EST C C SONEST ESTIMATES THE 1-NORM OF A SQUARE, REAL MATRIX A. C REVERSE COMMUNICATION IS USED FOR EVALUATING C MATRIX-VECTOR PRODUCTS. C C ON ENTRY C C N INTEGER C THE ORDER OF THE MATRIX. N .GE. 1. C C ISGN INTEGER(N) C USED AS WORKSPACE. C C KASE INTEGER C = 0. C C ON INTERMEDIATE RETURNS C C KASE = 1 OR 2. C C X REAL(N) C MUST BE OVERWRITTEN BY C C A*X, IF KASE=1, C TRANSPOSE(A)*X, IF KASE=2, C C AND SONEST MUST BE RE-CALLED, WITH ALL THE OTHER C PARAMETERS UNCHANGED. C C ON FINAL RETURN C C KASE = 0. C C EST REAL C CONTAINS AN ESTIMATE (A LOWER BOUND) FOR NORM(A). C C V REAL(N) C = A*W, WHERE EST = NORM(V)/NORM(W) C (W IS NOT RETURNED). C C THIS VERSION DATED MARCH 16, 1988. C NICK HIGHAM, UNIVERSITY OF MANCHESTER. C C REFERENCE C N.J. HIGHAM (1987) FORTRAN CODES FOR ESTIMATING C THE ONE-NORM OF A REAL OR COMPLEX MATRIX, WITH APPLICATIONS C TO CONDITION ESTIMATION, NUMERICAL ANALYSIS REPORT NO. 135, C UNIVERSITY OF MANCHESTER, MANCHESTER M13 9PL, ENGLAND. C C SUBROUTINES AND FUNCTIONS C BLAS ISAMAX, SASUM, SCOPY C GENERIC ABS, NINT, REAL, SIGN C PARAMETER (ITMAX = 5) PARAMETER (ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0) C C INTERNAL VARIABLES INTEGER I, ITER, J, JLAST, JUMP REAL ALTSGN, ESTOLD, TEMP C SAVE C IF (KASE .EQ. 0) THEN DO 10,I = 1,N X(I) = ONE/REAL(N) 10 CONTINUE KASE = 1 JUMP = 1 RETURN ENDIF C GOTO (100, 200, 300, 400, 500) JUMP C C ................ ENTRY (JUMP = 1) C FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. C 100 CONTINUE IF (N .EQ. 1) THEN V(1) = X(1) EST = ABS(V(1)) C ... QUIT GOTO 510 ENDIF EST = SASUM(N,X,1) C DO 110,I = 1,N X(I) = SIGN(ONE,X(I)) ISGN(I) = NINT(X(I)) 110 CONTINUE KASE = 2 JUMP = 2 RETURN C C ................ ENTRY (JUMP = 2) C FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. C 200 CONTINUE J = ISAMAX(N,X,1) ITER = 2 C C MAIN LOOP - ITERATIONS 2,3,...,ITMAX. C 220 CONTINUE DO 230,I = 1,N X(I) = ZERO 230 CONTINUE X(J) = ONE KASE = 1 JUMP = 3 RETURN C C ................ ENTRY (JUMP = 3) C X HAS BEEN OVERWRITTEN BY A*X. C 300 CONTINUE CALL SCOPY(N,X,1,V,1) ESTOLD = EST EST = SASUM(N,V,1) DO 310,I = 1,N IF ( NINT( SIGN(ONE,X(I)) ) .NE. ISGN(I) ) GOTO 320 310 CONTINUE C REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. GOTO 410 C 320 CONTINUE C TEST FOR CYCLING. IF (EST .LE. ESTOLD) GOTO 410 C DO 330,I = 1,N X(I) = SIGN(ONE,X(I)) ISGN(I) = NINT(X(I)) 330 CONTINUE KASE = 2 JUMP = 4 RETURN C C ................ ENTRY (JUMP = 4) C X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. C 400 CONTINUE JLAST = J J = ISAMAX(N,X,1) IF ( ( X(JLAST) .NE. ABS(X(J)) ) .AND. + (ITER .LT. ITMAX) ) THEN ITER = ITER + 1 GOTO 220 ENDIF C C ITERATION COMPLETE. FINAL STAGE. C 410 CONTINUE ALTSGN = ONE DO 420,I = 1,N X(I) = ALTSGN * (ONE + REAL(I-1)/REAL(N-1)) ALTSGN = -ALTSGN 420 CONTINUE KASE = 1 JUMP = 5 RETURN C C ................ ENTRY (JUMP = 5) C X HAS BEEN OVERWRITTEN BY A*X. C 500 CONTINUE TEMP = TWO*SASUM(N,X,1)/REAL(3*N) IF (TEMP. GT. EST) THEN CALL SCOPY(N,X,1,V,1) EST = TEMP ENDIF C 510 KASE = 0 RETURN C END PROGRAM SONRUN C C THIS IS A TEST DRIVER FOR SONEST. C IT APPLIES SONEST TO SEVERAL CLASSES OF MATRIX AND EVALUATES C VARIOUS PERFORMANCE MEASURES. IT SERVES THE FOLLOWING PURPOSES. C C (1) ILLUSTRATES HOW TO USE SONEST (WITH LINPACK'S SGEFA/SGESL). C C (2) RIGOROUSLY TESTS SONEST. EXERCISES ALL STATEMENTS AND C FEATURES OF THE CODE. C C (3) IS A TOOL FOR STUDYING THE BEHAVIOUR OF SONEST. C C RESULTS ARE REPORTED SEPARATELY FOR EACH OF FIVE CLASSES OF TEST C MATRIX. WITHIN EACH CLASS THE MINIMUM, MAXIMUM, AND AVERAGE VALUES C OF SEVEN STATISTICS ARE REPORTED FOR EACH COMBINATION OF CONDITION C NUMBER (WHERE APPLICABLE) AND N. C C TO CHECK THAT SONEST IS WORKING CORRECTLY CHECK THE FOLLOWING. C - STATISTIC 1 (SONEST ESTIMATE/TRUE NORM). C THE MINIMUM VALUES SHOULD GENERALLY BE AT LEAST 0.1. C THE POSSIBILITY OF SMALLER VALUES CANNOT BE RULED C OUT, SINCE THE ALGORITHM UNDERLYING THE CODE IS C NOT COMPLETELY RELIABLE. C - STATISTIC 4 (NULL VECTOR RATIO). C THE MINIMUM VALUE SHOULD BE 1.0, EXCEPT THAT THERE C WILL BE SOME DETERIORATION FOR LARGER CONDITION C NUMBERS, AS NOTED BELOW. C - THE SUMMARY OF THE NUMBER OF ITERATIONS. THIS SHOULD SHOW C RELATIVELY FEW OCCURRENCES OF 4 OR 5 ITERATIONS. C C THE OTHER STATISTICS ARE OF INTEREST FOR INVESTIGATING AND C UNDERSTANDING THE PERFORMANCE OF SONEST. FOR EXAMPLE, ONE C MIGHT EXAMINE THE FOLLOWING ASPECTS. C C - THE NUMBER OF TIMES THE EXTRA ESTIMATE IS USED. C - THE NUMBER OF TIMES A REPEATED LINEAR SYSTEM IS DETECTED ('SAVED C SOLVE'). C - THE FREQUENCY OF EXACT ESTIMATES (VALUE OF 1.0 FOR STATISTIC 1). C - HOW THE SONEST ESTIMATES COMPARE WITH THOSE OF SGECO (STATISTIC 6). C - HOW ALL THE ABOVE FACTORS VARY WITH N, THE CONDITION NUMBER, AND THE C TYPE OF MATRIX. C C THE USER SHOULD ADJUST THE PARAMETERS NCOND, NOFN, ITIMES, AND THE C ARRAY NVALS, TO TASTE - E.G. TO SUIT THE COMPUTER TIME AVAILABLE. C THE MODEST DEFAULT VALUES BELOW ARE SUITABLE FOR AN INITIAL TEST C RUN. FOR EFFICIENCY, THE CASES N=1, AND MODE=-2, SHOULD BE RUN C WITH ITIMES=1, SEPARATELY FROM THE MAIN TESTING WITH ITIMES>1 C C THIS VERSION DATED MARCH 16TH, 1988. C NICK HIGHAM, UNIVERSITY OF MANCHESTER. C C MACHINE DEPENDENCE C C THE ONLY MACHINE DEPENDENT PARAMETER IS EPS IN THE FOLLOWING C PARAMETER STATEMENT. SET THIS TO THE MACHINE PRECISION C (UNIT ROUNDOFF). PARAMETER (EPS=2.0**(-23)) C C SUBROUTINES AND FUNCTIONS. C C LINPACK - SGECO, SGEFA, SGESL, SGEDI. C BLAS - THOSE REQUIRED BY THE ABOVE LINPACK ROUTINES AND BY SONEST. C RANDOM NUMBER GENERATORS URAND AND RANDN, SUPPLIED WITH THIS CODE. C C THE TEST MATRICES FOR MODE.GE.0 WILL HAVE 2-NORM CONDITION NUMBERS C RANGING FROM 1 TO CONDMX, TAKING ON NCOND DIFFERENT VALUES THAT C ARE POWERS OF CONDPW (A REAL VARIABLE SET BELOW). C THE MAXIMUM CONDITION NUMBER, CONDMX, IS JUST BELOW THE LEVEL C CORRESPONDING TO THE MATRIX BEING 'SINGULAR TO WORKING PRECISION'. C C N.B. FOR CONDITION NUMBERS ON THE ORDER OF 1.0/EPS, ROUNDING ERRORS C ARE SIGNIFICANT AND THE RATIOS NEED TO BE INTERPRETED WITH CARE. C STATISTICS 4, 5 AND 6 (NULL VECTOR RATIO, SGECO/TRUE NORM C AND SGECO/SONEST) ARE ALL 'UNRELIABLE' FOR LARGE CONDITION C NUMBERS. HOWEVER, STATISTIC 1 (SONEST/TRUE NORM) DOES NOT C DETERIORATE FOR LARGE CONDITION NUMBERS BECAUSE OF CORRELATION C WITH THE 'TRUE' CONDITION NUMBER OBTAINED VIA SGEDI (ESSENTIALLY C THESE TWO ESTIMATES USE THE SAME COMPUTED COLUMNS OF INV(A)). C PARAMETER (CONDMX = (1.0/EPS)/10.0) PARAMETER (NCOND=3) C C NUMBER OF N VALUES PARAMETER (NOFN=3) C NUMBER OF RANDOM REPETITIONS FOR EACH COMBINATION OF MATRIX TYPE, C N, AND CONDITION NUMBER. PARAMETER (ITIMES=2) C C NUMBER OF STATISTICS TO INVESTIGATE PARAMETER (NESTS=7) C SONEST TAKES BETWEEN 1 AND ITMAX ITERATIONS, EXCEPT C IN THE TRIVIAL CASE N=1 WHEN IT TAKES 'HALF AN ITERATION' C (WHICH WE COUNT HERE AS ZERO ITERATIONS). PARAMETER (ITMAX=5) C C LDA MUST BE AT LEAST AS BIG AS THE MAXIMUM VALUE OF N. PARAMETER (LDA=50, LDASQ=LDA*LDA) C REAL A(LDA,LDA), X(LDA), V(LDA), ASAVE(LDA,LDA) REAL W1(LDA), W2(LDA), DET(2) REAL S(NESTS,ITIMES), RES(NESTS,NOFN,NCOND,3) INTEGER NVALS(NOFN), IPIVOT(LDA), ISGN(LDA), ITS(0:ITMAX) C C MULTIPLICATIVE FACTOR FOR GENERATING CONDITION NUMBERS. CONDPW = CONDMX**(1.0/REAL(NCOND-1)) C C THE NOFN DIFFERENT VALUES OF N TO BE USED. DATA NVALS/1,5,10/ * C INITIALISE A TO ZERO. DATA A/LDASQ*0.0/ C C FOR EACH OF THE FIVE CLASSES OF MATRIX... DO 60,MODE = -2,2 C ISING = 0 IF (MODE.LT.0) THEN C FOR MODE.LT.0 THE CONDITION NUMBER IS NOT A PARAMETER. INCOND = 1 ELSE INCOND = NCOND END IF C C RESET SONEST ITERATION STATISTICS. DO 1,I=0,ITMAX ITS(I) = 0 1 CONTINUE ISAVE = 0 IXTRA = 0 C C FOR EACH VALUE OF N... DO 30,I=1,NOFN N=NVALS(I) C FOR EACH CONDITION NUMBER... DO 30,J=1,INCOND C DO 10,K=1,ITIMES C SEED FOR THE RANDOM NUMBER GENERATOR. IRAND = K*J*I+1 3 IF (MODE.GE.0) THEN COND = CONDPW**(J-1) CALL UDV(N,A,LDA,COND,MODE,W1,W2,IRAND) ELSEIF (MODE.EQ.-1) THEN CALL RANDA (N,A,LDA,IRAND) ELSEIF (MODE.EQ.-2) THEN CALL BIDIAG (N,A,LDA) END IF CALL SCOPY(LDA*LDA,A,1,ASAVE,1) C ANORM = 0.0 DO 2,L=1,N ANORM = MAX(ANORM, SASUM(N,A(1,L),1)) 2 CONTINUE C GET PA=LU FACTORISATION TOGETHER WITH LINPACK CONDITION ESTIMATE. CALL SGECO(A,LDA,N,IPIVOT,RCOND,W1) C C IN THE UNLIKELY EVENT OF A SINGULAR U IN PA=LU (FOR MODE.GE.-1), C GENERATE ANOTHER A. IF (RCOND .EQ. 0) THEN ISING = ISING+1 GOTO 3 END IF C C ------------------------ C THIS SECTION OF CODE APPLIES SONEST. THE FIRST LINE, AND THE ONES C INVOLVING ISOLVE AND ESTOLD, COULD BE OMITTED IN PRACTICE; C THEY ARE USED HERE TO MONITOR THE PERFORMANCE OF SONEST. EST = 0.0 ISOLVE = 0 KASE = 0 5 CALL SONEST(N,V,X,ISGN,EST,KASE) IF (KASE .NE. 0) THEN ISOLVE = ISOLVE+1 ESTOLD = EST CALL SGESL(A,LDA,N,IPIVOT,X,KASE-1) GOTO 5 END IF C ------------------------ C C MAKE DEDUCTIONS ABOUT HOW MANY ITERATIONS (ITER) WERE USED, AND C WHETHER THERE WAS A 'SAVED SOLVE' DUE TO A REPEATED 'XI' VECTOR. C N.B. WE CAN'T DISTINGUISH BETWEEN DETECTION OF A REPEATED VECTOR C AND DETECTION OF CYCLING VIA THE EXPLICIT TEST FOR CYCLING. C HOWEVER, CYCLING IS EXTREMELY UNLIKELY (AND IT CANNOT HAPPEN IN C EXACT ARITHMETIC). MOREOVER, CYCLING CONSISTING OF MOVING BACK TO C THE MOST RECENT VERTEX WILL BE PICKED UP BY THE REPEATED VECTOR C TEST. IF (N .GT. 1) THEN IF (EST .NE. ESTOLD) IXTRA = IXTRA+1 ITSLVE = ISOLVE-1 ITER = INT( (ITSLVE+1)/2 ) IF ( MOD(ITSLVE,2) .EQ. 1 ) ISAVE = ISAVE+1 ITS(ITER) = ITS(ITER)+1 ELSE ITER = INT(ISOLVE/2) ITS(ITER) = ITS(ITER)+1 END IF C C CHECK APPROXIMATE NULL VECTOR. DO 6,L=1,N X(L) = SDOT(N,ASAVE(L,1),LDA,V,1) 6 CONTINUE XEST=SASUM(N,V,1)/SASUM(N,X,1) C C COMPUTE THE 'TRUE' VALUE OF NORM1(INV(A)). CALL SGEDI(A,LDA,N,IPIVOT,DET,W1,1) AINORM = 0.0 DO 7,L=1,N AINORM = MAX(AINORM,SASUM(N,A(1,L),1)) 7 CONTINUE COND1 = ANORM*AINORM C C FORM THE STATISTICS OF INTEREST. 9 S(1,K) = EST/AINORM S(2,K) = REAL(ITER) S(3,K) = REAL(ISOLVE) S(4,K) = XEST/EST S(5,K) = 1.0/(RCOND*COND1) S(6,K) = EST*ANORM*RCOND S(7,K) = COND1 C 10 CONTINUE C C RETAIN ONLY THE MINIMUM, MAXIMUM AND AVERAGE VALUES OF THE STATS. DO 20,L=1,NESTS SMIN=S(L,1) SMAX=S(L,1) SUM=0.0 DO 15,K=1,ITIMES SMIN=MIN(S(L,K),SMIN) SUM=SUM+S(L,K) SMAX=MAX(S(L,K),SMAX) 15 CONTINUE RES(L,I,J,1)=SMIN RES(L,I,J,2)=SMAX RES(L,I,J,3)=SUM/REAL(ITIMES) 20 CONTINUE 30 CONTINUE C C PRINT RESULTS WRITE(6,'(//)') WRITE(6,*) 'RESULTS FOR MODE = ',MODE,' ITIMES = ',ITIMES IF (MODE .EQ. -2) THEN WRITE(6,*) 'UPPER BIDIAGONAL, U(I,J) = 1' ELSE IF (MODE .EQ. -1) THEN WRITE(6,*) 'ELTS OF A FROM UNIFORM (-1,1) DISTRIBUTION' ELSE IF (MODE .EQ. 0) THEN WRITE(6,*) 'ONE LARGE SINGULAR VALUE' ELSE IF (MODE .EQ. 1) THEN WRITE(6,*) 'ONE SMALL SINGULAR VALUE' ELSE IF (MODE .EQ. 2) THEN WRITE(6,*) 'EXPONENTIALLY DISTRIBUTED SINGULAR VALUES' END IF WRITE(6,*) '##############################################' WRITE(6,*) WRITE(6,*) 'SONEST ITERATION BEHAVIOUR' WRITE(6,'(A18,6I6)') 'NO. ITERATIONS ', 0,1,2,3,4,5 WRITE(6,'(A18,6I6)') ' # TIMES', (ITS(I),I=0,ITMAX) WRITE(6,'(A36,I8)') ' NUMBER OF SAVED SOLVES (OR CYCLING)', ISAVE WRITE(6,'(A31,I8)') ' NUMBER OF TIMES EXTRA EST USED', IXTRA WRITE(6,'(A25,I8)') ' TOTAL NUMBER OF MATRICES ', + NOFN*INCOND*ITIMES WRITE(6,'(A22,I8/)') ' NUMBER OF SINGULAR U =', ISING DO 50,L=1,NESTS WRITE(6,*) WRITE(6,*) '********************************************' WRITE(6,*) 'STATISTIC NO ',L IF (L .EQ. 1) THEN WRITE(6,*) 'RATIO - SONEST ESTIMATE / NORM1(INV(A))' ELSE IF (L .EQ. 2) THEN WRITE(6,*) 'NUMBER OF ITERATIONS' ELSE IF (L .EQ. 3) THEN WRITE(6,*) 'NUMBER OF LINEAR SYSTEMS SOLVED' ELSE IF (L .EQ. 4) THEN WRITE(6,*) 'CHECK ON APPROX NULL VECTOR - SHOULD BE 1.0' ELSE IF (L .EQ. 5) THEN WRITE(6,*) 'RATIO - LINPACK ESTIMATE / NORM1(INV(A))' ELSE IF (L .EQ. 6) THEN WRITE(6,*) 'RATIO - SONEST ESTIMATE / LINPACK ESTIMATE' ELSE IF (L .EQ. 7) THEN WRITE(6,*) '1-NORM CONDITION NUMBER' END IF WRITE(6,*) 'MIN, MAX THEN AVE' DO 40,K=1,3 WRITE(6,*) IF (MODE .LE. -1) THEN WRITE(6,'(A5,6I12)') ' N ', (NVALS(I),I=1,NOFN) ELSE WRITE(6,'(A12,6I12)') 'COND / N ',(NVALS(I),I=1,NOFN) END IF WRITE(6,*) '--------------------------------------' DO 35,J=1,INCOND IF (MODE .LE. -1) THEN WRITE(6,'(5X,6E12.4)') (RES(L,I,J,K),I=1,NOFN) ELSE WRITE(6,'(7E12.4)') CONDPW**(J-1),(RES(L,I,J,K),I=1,NOFN) END IF 35 CONTINUE 40 CONTINUE 50 CONTINUE 60 CONTINUE C END C SUBROUTINE UDV(N,A,LDA,COND2,MODE,X,D,IRAND) C GENERATES RANDOM A WITH PRE-ASSIGNED SINGULAR VALUE DISTRIBUTION. C A = UDV, WHERE U AND V ARE RANDOM ORTHOGONAL MATRICES AND C D IS DIAGONAL WITH THE SINGULAR VALUES ON ITS DIAGONAL. C REAL A(LDA,N), X(N), D(N) DO 5,I=1,N DO 5,J=1,N A(I,J) = 0.0 5 CONTINUE IF (MODE.EQ.0) THEN A(1,1) = 1.0 DO 10,I=2,N A(I,I) = 1.0/COND2 10 CONTINUE ELSE IF (MODE.EQ.1) THEN DO 15,I=1,N-1 A(I,I) = 1.0 15 CONTINUE A(N,N) = 1.0/COND2 ELSE IF (MODE.EQ.2) THEN A(1,1) = 1.0 IF (N.GT.1) THEN ALPHA = COND2**(-1.0/REAL(N-1)) DO 20,I=2,N A(I,I) = ALPHA*A(I-1,I-1) 20 CONTINUE END IF END IF CALL UTIMEA(N,A,LDA,X,D,1,IRAND) C TRANSPOSE A C (AVOIDS NEED FOR AN ANALOGUE OF UTIMEA FOR POST-MULTIPLICATION). DO 30,I=1,N-1 DO 30,J=I+1,N T = A(I,J) A(I,J) = A(J,I) A(J,I) = T 30 CONTINUE CALL UTIMEA(N,A,LDA,X,D,2,IRAND) END C SUBROUTINE UTIMEA(N,A,LDA,X,D,ICALL,IRAND) C PRE-MULTIPLIES A BY A RANDOM ORTHOGONAL MATRIX U, C OVERWRITING A. U IS GENERATED USING THE METHOD OF G.W. STEWART C (SIAM J. NUMER. ANAL. 17, 1980, PP. 403-409). C REAL A(LDA,N), X(N), D(N), RANDN INTEGER COL, ROW DO 40,I=2,N K = N-I+1 S = 0.0 DO 10,J=K,N T = RANDN(IRAND) X(J) = T S = S+T*T 10 CONTINUE S = SQRT(S) S = SIGN(S,X(K)) D(K) = SIGN(1.0,-S) PI = S*(S+X(K)) X(K) = X(K)+S IF (ICALL.EQ.1) THEN J = K ELSE J = 1 END IF DO 35,COL=J,N SP = 0.0 DO 20,ROW=K,N SP = SP+X(ROW)*A(ROW,COL) 20 CONTINUE DO 30,ROW=K,N A(ROW,COL) = A(ROW,COL)-SP*X(ROW)/PI 30 CONTINUE 35 CONTINUE 40 CONTINUE D(N) = SIGN(1.0,RANDN(IRAND)) DO 50,I=1,N T = D(I) DO 50,J=1,N A(I,J) = T*A(I,J) 50 CONTINUE END C SUBROUTINE RANDA (N,A,LDA,IRAND) C GENERATES RANDOM A WITH ELEMENTS FROM UNIFORM (-1,1) DISTRIBUTION. REAL A(LDA,N),URAND DO 10,J=1,N DO 10,I=1,N A(I,J) = -1.0 + 2.0*URAND(IRAND) 10 CONTINUE END C SUBROUTINE BIDIAG (N,A,LDA) REAL A(LDA,N) DO 10,J=1,N DO 10,I=1,N IF ( (J.EQ.I) .OR. (J.EQ.I+1) ) THEN A(I,J) = 1.0 ELSE A(I,J) = 0.0 END IF 10 CONTINUE END SUBROUTINE CONEST (N, V, X, EST, KASE) INTEGER N, KASE COMPLEX V(N), X(N) REAL EST C C CONEST ESTIMATES THE 1-NORM OF A SQUARE, COMPLEX MATRIX A. C REVERSE COMMUNICATION IS USED FOR EVALUATING C MATRIX-VECTOR PRODUCTS. C C ON ENTRY C C N INTEGER C THE ORDER OF THE MATRIX. N .GE. 1. C C KASE INTEGER C = 0. C C ON INTERMEDIATE RETURNS C C KASE = 1 OR 2. C C X COMPLEX(N) C MUST BE OVERWRITTEN BY C C A*X, IF KASE=1, C CTRANS(A)*X, IF KASE=2, C C WHERE CTRANS(A) IS THE CONJUGATE TRANSPOSE, C AND CONEST MUST BE RE-CALLED, WITH ALL THE OTHER C PARAMETERS UNCHANGED. C C ON FINAL RETURN C C KASE = 0. C C EST REAL C CONTAINS AN ESTIMATE (A LOWER BOUND) FOR NORM(A). C C V COMPLEX(N) C = A*W, WHERE EST = NORM(V)/NORM(W) C (W IS NOT RETURNED). C C THIS VERSION DATED MARCH 16, 1988. C NICK HIGHAM, UNIVERSITY OF MANCHESTER. C C REFERENCE C N.J. HIGHAM (1987) FORTRAN CODES FOR ESTIMATING C THE ONE-NORM OF A REAL OR COMPLEX MATRIX, WITH APPLICATIONS C TO CONDITION ESTIMATION, NUMERICAL ANALYSIS REPORT NO. 135, C UNIVERSITY OF MANCHESTER, MANCHESTER M13 9PL, ENGLAND. C C SUBROUTINES AND FUNCTIONS C BLAS CCOPY C MODIFIED BLAS ICMAX1, SCSUM1 C GENERIC ABS, CMPLX, REAL C PARAMETER (ITMAX = 5) PARAMETER (ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0) COMPLEX CZERO, CONE PARAMETER (CZERO = (0.0E0, 0.0E0), CONE = (1.0E0, 0.0E0)) * C C INTERNAL VARIABLES INTEGER I, ITER, J, JLAST, JUMP REAL ALTSGN, ESTOLD, TEMP C SAVE C IF (KASE .EQ. 0) THEN DO 10,I = 1,N X(I) = CMPLX( ONE/REAL(N) ) 10 CONTINUE KASE = 1 JUMP = 1 RETURN ENDIF C GOTO (100, 200, 300, 400, 500) JUMP C C ................ ENTRY (JUMP = 1) C FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. C 100 CONTINUE IF (N .EQ. 1) THEN V(1) = X(1) EST = ABS(V(1)) C ... QUIT GOTO 510 ENDIF EST = SCSUM1(N,X,1) C DO 110,I = 1,N IF ( ABS(X(I)) .NE. ZERO ) THEN X(I) = X(I) / CMPLX( ABS(X(I)) ) ELSE X(I) = CONE ENDIF 110 CONTINUE KASE = 2 JUMP = 2 RETURN C C ................ ENTRY (JUMP = 2) C FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X. C 200 CONTINUE J = ICMAX1(N,X,1) ITER = 2 C C MAIN LOOP - ITERATIONS 2,3,...,ITMAX. C 220 CONTINUE DO 230,I = 1,N X(I) = CZERO 230 CONTINUE X(J) = CONE KASE = 1 JUMP = 3 RETURN C C ................ ENTRY (JUMP = 3) C X HAS BEEN OVERWRITTEN BY A*X. C 300 CONTINUE CALL CCOPY(N,X,1,V,1) ESTOLD = EST EST = SCSUM1(N,V,1) C C TEST FOR CYCLING. IF (EST .LE. ESTOLD) GOTO 410 C DO 330,I = 1,N IF ( ABS(X(I)) .NE. ZERO ) THEN X(I) = X(I) / CMPLX( ABS(X(I)) ) ELSE X(I) = CONE ENDIF 330 CONTINUE KASE = 2 JUMP = 4 RETURN C C ................ ENTRY (JUMP = 4) C X HAS BEEN OVERWRITTEN BY CTRANS(A)*X. C 400 CONTINUE JLAST = J J = ICMAX1(N,X,1) IF ( ( REAL(X(JLAST)) .NE. ABS( REAL(X(J)) ) ) .AND. + (ITER .LT. ITMAX) ) THEN ITER = ITER + 1 GOTO 220 ENDIF C C ITERATION COMPLETE. FINAL STAGE. C 410 CONTINUE ALTSGN = ONE DO 420,I = 1,N X(I) = CMPLX( ALTSGN * (ONE + REAL(I-1)/REAL(N-1)) ) ALTSGN = -ALTSGN 420 CONTINUE KASE = 1 JUMP = 5 RETURN C C ................ ENTRY (JUMP = 5) C X HAS BEEN OVERWRITTEN BY A*X. C 500 CONTINUE TEMP = TWO*SCSUM1(N,X,1)/REAL(3*N) IF (TEMP. GT. EST) THEN CALL CCOPY(N,X,1,V,1) EST = TEMP ENDIF C 510 KASE = 0 RETURN C END