program cyldr real kl,ka,er,ei,theta real kdel,pi,degrees real kz(200),kr,ku,kzn,kznn real reca,rect real maga,magt,phsa,phst integer n,dima,dimt,orda,ordt integer i,j,piva(200),pivt(200) complex e,ii complex eia(200),eit(200),za(200,200),zt(200,200) complex wrka(200),wrkt(200) read *, kl read *, ka read *, er read *, ei read *, theta read *, n kdel=kl/(2.0*n) ii=cmplx(0.0,1.0) pi=4.0*atan(1.0) degrees=pi/180. e=cmplx(er,ei) dima=200 dimt=200 orda=n ordt=n open(unit=23,file='cyldr.grg') print *, 'eia(i)' do 5 i=1,n kz(i)=kl*(2.0*i-n-1)/(2.0*n) eia(i)=cexp(ii*cos(theta*degrees)*kz(i)) eit(i)=eia(i) kzn=kz(i)/kl kznn=(0.5*kl+kz(i))/ka write(23,*) kz(i),kzn,kznn,cabs(eia(i)), % real(eia(i)),aimag(eia(i)) print *, i,eia(i) 5 continue print*, 'j,za(i,j),zt(i,j)*(e+1.0)/2.0,kr,kdel' do 10 i=1,n do 20 j=1,n if (i .eq. j) then kr=sqrt(ka**2+kdel**2) za(i,j)=1.0-0.5*(e-1.0)*(ii*ka**2*kdel+kdel*(kr-kdel) & +(ka**2/2.0)*alog(kr+kdel)-(ka**2/2.0)*alog(kr-kdel)- & cexp(ii*kdel)*2.0 + cexp(ii*kr)*2.0*kdel/kr) c zt(i,j)=(2.0/(e+1.0)) c & -(e-1.0)/(e+1.0)*(ii*2.0*ka**2*kdel/3.0-kdel**3 c & -kdel/kr-(ka**2/4.0)*alog(kr-kdel) c & +(ka**2/4.0)*alog(kr+kdel)) zt(i,j)=(2.0/(e+1.0)) & -(e-1.0)/(e+1.0)*(-kdel**2 + & kdel*kr - kdel/kr + ii*2.0*ka**2*kdel/3.0) else ku=(kz(i)-kz(j)) kr=sqrt(ka**2+ku**2) za(i,j)=(e-1.0)*(-0.5)*(ii*2.0*kdel*(cexp(ii*abs(ku))- & cexp(ii*kr)) + kdel**3*ka**2*cexp(ii*kr)/(3.0*kr**3) & -ii*kdel**3*ka**2*cexp(ii*kr)/(3.0*kr**2) & +cexp(ii*abs(ku-kdel))*(ku-kdel)/abs(ku-kdel) & -cexp(ii*sqrt(ka**2+(ku-kdel)**2))*(ku-kdel)/ & sqrt(ka**2+(ku-kdel)**2) & -cexp(ii*abs(ku+kdel))*(ku+kdel)/abs(ku+kdel) & +cexp(ii*sqrt(ka**2+(ku+kdel)**2))*(ku+kdel)/ & sqrt(ka**2+(ku+kdel)**2)) c zt(i,j)=(1.0-e)/(e+1.0)*(2.0*kdel*(ii*cexp(ii*abs(ku)) c & -ii*cexp(ii*kr)+(ii*ka**2*0.5)*cexp(ii*kr)/kr**2 c & -(ka**2*0.5)*cexp(ii*kr)/kr**3) c & +kdel**3*(2.0*ii*ka**2*ku**2*cexp(ii*kr)/kr**6 c & -2.0*ka**2*ku**2*cexp(ii*kr)/kr**7)) zt(i,j)=(1.0-e)/(e+1.0)*(ii*2.0*kdel*(cexp(ii*abs(ku)) & -cexp(ii*kr)) + ka**2*kdel**3*cexp(ii*kr)/(3.0*kr**3) & -ii*ka**2*kdel**3*cexp(ii*kr)/(3.0*kr**2) & +0.5*ii*ka*cexp(ii*kr)*(atan((kdel-ku)/ka)+ & atan((kdel+ku)/ka)) -0.5*cexp(ii*kr)* & ((kdel-ku)/sqrt(ka**2+(kdel-ku)**2)+ & (kdel+ku)/sqrt(ka**2+(kdel+ku)**2))) endif if(i .eq. 75)then print*, j,za(i,j),zt(i,j)*(e+1.0)/2.0,kr,kdel else continue endif 20 continue 10 continue do 25 i=1,n c print *, i,za(i,50),zt(i,50) 25 continue c print *, 'ok1' call cgeco(za,dima,orda,piva,reca,wrka) c print *, 'ok2' call cgesl(za,dima,orda,piva,eia,0) c print *, 'ok3' call cgeco(zt,dimt,ordt,pivt,rect,wrkt) c print *, 'ok4' call cgesl(zt,dimt,ordt,pivt,eit,0) c print *, 'ok5' open(unit=20,file='cyldra.out') open(unit=21,file='cyldrt.out') pa=0.0 pt=0.0 do 30 i=1,n c print *, eia(i) maga=cabs(eia(i)) phsa=atan2(real(eia(i)),aimag(eia(i))) magt=cabs(eit(i)) phst=atan2(real(eit(i)),aimag(eit(i))) kzn=kz(i)/kl kznn=(0.5*kl+kz(i))/ka write(20,*)kz(i),kzn,kznn,maga,real(eia(i)),aimag(eia(i)) write(21,*)kz(i),kzn,kznn,magt,real(eit(i)),aimag(eit(i)) pa=pa+maga/n pt=pt+magt/n 30 continue write(20,*) pa write(21,*) pt print*, 'pa = ',pa,' pt= ',pt close(20) close(21) close(23) stop end 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 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 CCC 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 CCC 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 CCC 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 CCC 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 CCC 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 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 CCC 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 CCC Zero pivot implies this column already triangularized C IF (CABS1(A(L,K)) .EQ. 0.0E0) GO TO 40 C CCC 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 CCC Compute multipliers C T = -CMPLX(1.0E0,0.0E0)/A(K,K) CALL CSCAL(N-K,T,A(K+1,K),1) C CCC 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 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 CCC Code for unequal increments or equal increments CCC 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 CCC 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 CCC Code for unequal increments or equal increments CCC 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 CCC 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 CCC 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 CCC 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 CCC 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 CCC 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 CCC 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 CCC 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 CCC 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 CCC 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