C******************************************************************C C C C The following subroutines are standard LINPACK routines C C to perform L-U decomposition and back substitution on a C C single precision complex matrix. See CC-Memo 407 sec 2.1 C C for documentation on these routines. C C C C******************************************************************C C C SUBROUTINE CGECO(A,LDA,N,IPVT,RCOND,Z) C C C******************************************************************C C C C NAASA 2.1.042 CGECO FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C INTEGER LDA,N,IPVT(1) COMPLEX A(LDA,1),Z(1) 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 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 07/14/77 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 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 = CMPLX(1.0E0,0.0E0) DO 20 J = 1, N Z(J) = CMPLX(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 = CMPLX(1.0E0,0.0E0) WKM = CMPLX(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 = V 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) = CMPLX(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 C C C******************************************************************C C C SUBROUTINE CGEFA(A,LDA,N,IPVT,INFO) C C C******************************************************************C C C C NAASA 2.1.043 CGEFA FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C INTEGER LDA,N,IPVT(1),INFO COMPLEX A(LDA,1) 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 07/14/77 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C BLAS CAXPY,CSCAL,ICAMAX C FORTRAN ABS,AIMAG,CMPLX,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 = -CMPLX(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 C C C******************************************************************C C C SUBROUTINE CGESL(A,LDA,N,IPVT,B,JOB) C C C******************************************************************C C C C NAASA 2.1.044 CGESL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C INTEGER LDA,N,IPVT(1),JOB COMPLEX A(LDA,1),B(1) 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 07/14/77 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 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 C C NAASA 2.1.045 CGEDI FTN-A 05-02-78 THE UNIV OF MICH COMP CTR SUBROUTINE CGEDI(A,LDA,N,IPVT,DET,WORK,JOB) INTEGER LDA,N,IPVT(1),JOB COMPLEX A(LDA,1),DET(1),WORK(1) 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 07/14/77 . C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 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) = CMPLX(1.0E0,0.0E0) DET(2) = CMPLX(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) - CMPLX(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) + CMPLX(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) = CMPLX(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) = CMPLX(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) = CMPLX(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 C C******************************************************************C C C SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY) C C C******************************************************************C C C C NAASA 1.1.014 CAXPY FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1),CY(1),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)GOTO 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 C C C******************************************************************C COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY) C******************************************************************C C C C NAASA 1.1.012 CDOTC FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C FORMS THE DOT PRODUCT OF TWO VECTORS, CONJUGATING THE FIRST C VECTOR. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1),CY(1),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)GOTO 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 C C******************************************************************C C C SUBROUTINE CSCAL(N,CA,CX,INCX) C C C******************************************************************C C C C NAASA 1.1.019 CSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C SCALES A VECTOR BY A CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CA,CX(1) INTEGER I,INCX,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 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 C C C******************************************************************C C C SUBROUTINE CSSCAL(N,SA,CX,INCX) C C C******************************************************************C C C C NAASA 1.1.018 CSSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C SCALES A COMPLEX VECTOR BY A REAL CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL SA INTEGER I,INCX,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 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 C C C******************************************************************C INTEGER FUNCTION ICAMAX(N,CX,INCX) C******************************************************************C C C C NAASA 1.1.021 ICAMAX FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL SMAX INTEGER I,INCX,IX,N COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C ICAMAX = 1 IF(N.LE.1)RETURN IF(INCX.EQ.1)GOTO 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 C C C******************************************************************C REAL FUNCTION SCASUM(N,CX,INCX) C******************************************************************C C C C NAASA 1.1.010 SCASUM FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND C RETURNS A SINGLE PRECISION RESULT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL STEMP INTEGER I,INCX,N,NINCX C SCASUM = 0.0E0 STEMP = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GOTO 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 CCOPY(N,CX,INCX,CY,INCY) 00000010 C 00000020 C COPIES A VECTOR, X, TO A VECTOR, Y. 00000030 C JACK DONGARRA, LINPACK, 3/11/78. 00000040 C 00000050 COMPLEX CX(1),CY(1) 00000060 INTEGER I,INCX,INCY,IX,IY,N 00000070 C 00000080 IF(N.LE.0)RETURN 00000090 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 00000100 C 00000110 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 00000120 C NOT EQUAL TO 1 00000130 C 00000140 IX = 1 00000150 IY = 1 00000160 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 00000170 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 00000180 DO 10 I = 1,N 00000190 CY(IY) = CX(IX) 00000200 IX = IX + INCX 00000210 IY = IY + INCY 00000220 10 CONTINUE 00000230 RETURN 00000240 C 00000250 C CODE FOR BOTH INCREMENTS EQUAL TO 1 00000260 C 00000270 20 DO 30 I = 1,N 00000280 CY(I) = CX(I) 00000290 30 CONTINUE 00000300 RETURN 00000310 END 00000320 SUBROUTINE CROTG(CA,CB,C,S) 00000010 COMPLEX CA,CB,S 00000020 REAL C 00000030 REAL NORM,SCALE 00000040 COMPLEX ALPHA 00000050 IF (CABS(CA) .NE. 0.) GO TO 10 00000060 C = 0. 00000070 S = (1.,0.) 00000080 CA = CB 00000090 GO TO 20 00000100 10 CONTINUE 00000110 SCALE = CABS(CA) + CABS(CB) 00000120 NORM = SCALE * SQRT((CABS(CA/SCALE))**2 + (CABS(CB/SCALE))**2) 00000130 ALPHA = CA /CABS(CA) 00000140 C = CABS(CA) / NORM 00000150 S = ALPHA * CONJG(CB) / NORM 00000160 CA = ALPHA * NORM 00000170 20 CONTINUE 00000180 RETURN 00000190 END 00000200 SUBROUTINE CSROT (N,CX,INCX,CY,INCY,C,S) 00000010 C 00000020 C APPLIES A PLANE ROTATION, WHERE THE COS AND SIN (C AND S) ARE REAL00000030 C AND THE VECTORS CX AND CY ARE COMPLEX. 00000040 C JACK DONGARRA, LINPACK, 3/11/78. 00000050 C 00000060 COMPLEX CX(1),CY(1),CTEMP 00000070 REAL C,S 00000080 INTEGER I,INCX,INCY,IX,IY,N 00000090 C 00000100 IF(N.LE.0)RETURN 00000110 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 00000120 C 00000130 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL 00000140 C TO 1 00000150 C 00000160 IX = 1 00000170 IY = 1 00000180 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 00000190 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 00000200 DO 10 I = 1,N 00000210 CTEMP = C*CX(IX) + S*CY(IY) 00000220 CY(IY) = C*CY(IY) - S*CX(IX) 00000230 CX(IX) = CTEMP 00000240 IX = IX + INCX 00000250 IY = IY + INCY 00000260 10 CONTINUE 00000270 RETURN 00000280 C 00000290 C CODE FOR BOTH INCREMENTS EQUAL TO 1 00000300 C 00000310 20 DO 30 I = 1,N 00000320 CTEMP = C*CX(I) + S*CY(I) 00000330 CY(I) = C*CY(I) - S*CX(I) 00000340 CX(I) = CTEMP 00000350 30 CONTINUE 00000360 RETURN 00000370 END 00000380 SUBROUTINE CSWAP (N,CX,INCX,CY,INCY) 00000010 C 00000020 C INTERCHANGES TWO VECTORS. 00000030 C JACK DONGARRA, LINPACK, 3/11/78. 00000040 C 00000050 COMPLEX CX(1),CY(1),CTEMP 00000060 INTEGER I,INCX,INCY,IX,IY,N 00000070 C 00000080 IF(N.LE.0)RETURN 00000090 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 00000100 C 00000110 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL 00000120 C TO 1 00000130 C 00000140 IX = 1 00000150 IY = 1 00000160 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 00000170 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 00000180 DO 10 I = 1,N 00000190 CTEMP = CX(IX) 00000200 CX(IX) = CY(IY) 00000210 CY(IY) = CTEMP 00000220 IX = IX + INCX 00000230 IY = IY + INCY 00000240 10 CONTINUE 00000250 RETURN 00000260 C 00000270 C CODE FOR BOTH INCREMENTS EQUAL TO 1 00000280 20 DO 30 I = 1,N 00000290 CTEMP = CX(I) 00000300 CX(I) = CY(I) 00000310 CY(I) = CTEMP 00000320 30 CONTINUE 00000330 RETURN 00000340 END 00000350 REAL FUNCTION SCNRM2( N, CX, INCX) 00000010 LOGICAL IMAG, SCALE 00000020 INTEGER NEXT 00000030 REAL CUTLO, CUTHI, HITEST, SUM, XMAX, ABSX, ZERO, ONE 00000040 COMPLEX CX(1) 00000050 DATA ZERO, ONE /0.0E0, 1.0E0/ 00000060 C 00000070 C UNITARY NORM OF THE COMPLEX N-VECTOR STORED IN CX() WITH STORAGE 00000080 C INCREMENT INCX . 00000090 C IF N .LE. 0 RETURN WITH RESULT = 0. 00000100 C IF N .GE. 1 THEN INCX MUST BE .GE. 1 00000110 C 00000120 C C.L.LAWSON , 1978 JAN 08 00000130 C 00000140 C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE 00000150 C HOPEFULLY APPLICABLE TO ALL MACHINES. 00000160 C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES. 00000170 C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES. 00000180 C WHERE 00000190 C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. 00000200 C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) 00000210 C V = LARGEST NO. (OVERFLOW LIMIT) 00000220 C 00000230 C BRIEF OUTLINE OF ALGORITHM.. 00000240 C 00000250 C PHASE 1 SCANS ZERO COMPONENTS. 00000260 C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO 00000270 C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO 00000280 C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M 00000290 C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. 00000300 C 00000310 C VALUES FOR CUTLO AND CUTHI.. 00000320 C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER 00000330 C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. 00000340 C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE00000350 C UNIVAC AND DEC AT 2**(-103) 00000360 C THUS CUTLO = 2**(-51) = 4.44089E-16 00000370 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. 00000380 C THUS CUTHI = 2**(63.5) = 1.30438E19 00000390 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. 00000400 C THUS CUTLO = 2**(-33.5) = 8.23181D-11 00000410 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 00000420 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / 00000430 C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / 00000440 DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / 00000450 C 00000460 IF(N .GT. 0) GO TO 10 00000470 SCNRM2 = ZERO 00000480 GO TO 300 00000490 C 00000500 10 ASSIGN 30 TO NEXT 00000510 SUM = ZERO 00000520 NN = N * INCX 00000530 C BEGIN MAIN LOOP 00000540 DO 210 I=1,NN,INCX 00000550 ABSX = ABS(REAL(CX(I))) 00000560 IMAG = .FALSE. 00000570 GO TO NEXT,(30, 50, 70, 90, 110) 00000580 30 IF( ABSX .GT. CUTLO) GO TO 85 00000590 ASSIGN 50 TO NEXT 00000600 SCALE = .FALSE. 00000610 C 00000620 C PHASE 1. SUM IS ZERO 00000630 C 00000640 50 IF( ABSX .EQ. ZERO) GO TO 200 00000650 IF( ABSX .GT. CUTLO) GO TO 85 00000660 C 00000670 C PREPARE FOR PHASE 2. 00000680 ASSIGN 70 TO NEXT 00000690 GO TO 105 00000700 C 00000710 C PREPARE FOR PHASE 4. 00000720 C 00000730 100 ASSIGN 110 TO NEXT 00000740 SUM = (SUM / ABSX) / ABSX 00000750 105 SCALE = .TRUE. 00000760 XMAX = ABSX 00000770 GO TO 115 00000780 C 00000790 C PHASE 2. SUM IS SMALL. 00000800 C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. 00000810 C 00000820 70 IF( ABSX .GT. CUTLO ) GO TO 75 00000830 C 00000840 C COMMON CODE FOR PHASES 2 AND 4. 00000850 C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.00000860 C 00000870 110 IF( ABSX .LE. XMAX ) GO TO 115 00000880 SUM = ONE + SUM * (XMAX / ABSX)**2 00000890 XMAX = ABSX 00000900 GO TO 200 00000910 C 00000920 115 SUM = SUM + (ABSX/XMAX)**2 00000930 GO TO 200 00000940 C 00000950 C 00000960 C PREPARE FOR PHASE 3. 00000970 C 00000980 75 SUM = (SUM * XMAX) * XMAX 00000990 C 00001000 85 ASSIGN 90 TO NEXT 00001010 SCALE = .FALSE. 00001020 C 00001030 C FOR REAL OR D.P. SET HITEST = CUTHI/N 00001040 C FOR COMPLEX SET HITEST = CUTHI/(2*N) 00001050 C 00001060 HITEST = CUTHI/FLOAT( N ) 00001070 C 00001080 C PHASE 3. SUM IS MID-RANGE. NO SCALING. 00001090 C 00001100 90 IF(ABSX .GE. HITEST) GO TO 100 00001110 SUM = SUM + ABSX**2 00001120 200 CONTINUE 00001130 C CONTROL SELECTION OF REAL AND IMAGINARY PARTS. 00001140 C 00001150 IF(IMAG) GO TO 210 00001160 ABSX = ABS(AIMAG(CX(I))) 00001170 IMAG = .TRUE. 00001180 GO TO NEXT,( 50, 70, 90, 110 ) 00001190 C 00001200 210 CONTINUE 00001210 C 00001220 C END OF MAIN LOOP. 00001230 C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. 00001240 C 00001250 SCNRM2 = SQRT(SUM) 00001260 IF(SCALE) SCNRM2 = SCNRM2 * XMAX 00001270 300 CONTINUE 00001280 RETURN 00001290 END 00001300 SUBROUTINE SROTG(SA,SB,C,S) 00000010 C 00000020 C CONSTRUCT GIVENS PLANE ROTATION. 00000030 C JACK DONGARRA, LINPACK, 3/11/78. 00000040 C 00000050 REAL SA,SB,C,S,ROE,SCALE,R,Z 00000060 C 00000070 ROE = SB 00000080 IF( ABS(SA) .GT. ABS(SB) ) ROE = SA 00000090 SCALE = ABS(SA) + ABS(SB) 00000100 IF( SCALE .NE. 0.0 ) GO TO 10 00000110 C = 1.0 00000120 S = 0.0 00000130 R = 0.0 00000140 GO TO 20 00000150 10 R = SCALE*SQRT((SA/SCALE)**2 + (SB/SCALE)**2) 00000160 R = SIGN(1.0,ROE)*R 00000170 C = SA/R 00000180 S = SB/R 00000190 20 Z = 1.0 00000200 IF( ABS(SA) .GT. ABS(SB) ) Z = S 00000210 IF( ABS(SB) .GE. ABS(SA) .AND. C .NE. 0.0 ) Z = 1.0/C 00000220 SA = R 00000230 SB = Z 00000240 RETURN 00000250 END 00000260