program circ real karc,ka,er,ei,theta,phi real kdel,pi,degrees real kphi(600),kr,ku,kl real reca,dd real eaha,eahb,eahc real eava,eavb,eavc integer n,dima,dimt,orda,ordt integer i,j,piv(600),pivt(600) complex e,ii complex phase(200),a,b,c,d,g complex zxx,zxy,zyy,zyx,zzz complex z(600,600),zt(600,600) complex eihx,eihy,eihz complex eivx,eivy,eivz complex eih(600),eiv(600) complex eah,ebh,ech complex eav,ebv,ecv complex wrka(600) read *, karc read *, kr read *, ka read *, er read *, ei read *, theta read *, phi kdel=ka/kr n=int(karc/(2.0*kdel)) kdel=karc/(2.0*n) ii=cmplx(0.0,1.0) pi=4.0*atan(1.0) degrees=pi/180. e=cmplx(er,ei) dima=600 dimt=600 orda=3*n ordt=3*n if (orda .gt. dima) then print*, 'WARNING, # of elements exceeds dimension' goto 100 else continue endif print*, 'n,kdel,ka' print*, n,kdel,ka do 5 i=1,n kphi(i)=karc*(2.0*i-n-1)/(2.0*n) phase(i)=cexp(ii*kr*sin(theta*degrees)*cos(kphi(i)-phi*degrees)) c eihx= -sin(phi*degrees)*phase(i) eihx=cmplx(0.0,0.0) eihy= cos(phi*degrees)*phase(i) eihz=cmplx(0.0,0.0) c eivx= cos(phi*degrees)*cos(theta*degrees)*phase(i) c eivy= sin(phi*degrees)*cos(theta*degrees)*phase(i) c eivz= -sin(theta*degrees)*phase(i) eivx=cmplx(0.0,0.0) eivy=cmplx(0.0,0.0) eivz=cmplx(0.0,0.0) eih(3*i-2)=eihx eih(3*i-1)=eihy eih(3*i)=eihz eiv(3*i-2)=eivx eiv(3*i-1)=eivy eiv(3*i)=eivz if(i .eq. 75)then print*, 'eihx,eihy,eihz,eivx,eivy,eivz' print*, eihx,eihy,eihz,eivx,eivy,eivz else continue endif 5 continue print*, 'j,ku,zyy,zxy,g,kr,kdel,kl' do 10 i=1,n do 20 j=1,n if (i .eq. j) then kl=sqrt(ka**2+(kdel*kr)**2) g=1.0-0.5*(e-1.0)*(ii*ka**2*(kdel*kr)+(kdel*kr) & *(kl-(kdel*kr)) & +(ka**2/2.0)*alog(kl+(kdel*kr))-(ka**2/2.0)* & alog(kl-(kdel*kr))- & cexp(ii*(kdel*kr))*2.0 + cexp(ii*kl)*2.0*(kdel*kr)/kl) d=1.0 & -((e-1.0)/2.0)*(-(kdel*kr)**2 + & (kdel*kr)*kl - (kdel*kr)/kl + ii*2.0*ka**2*(kdel*kr)/3.0) c & -((e-1.0)/2.0)*(ii*2.0*ka**2*(kdel*kr)/3.0-(kdel*kr)**3 c & -(kdel*kr)/kl-(ka**2/4.0)*alog(kl-(kdel*kr)) c & +(ka**2/4.0)*alog(kl+(kdel*kr))) zxx=d*cos(kphi(j))**2+g*sin(kphi(j))**2 zxy=cos(kphi(j))*sin(kphi(j))*(d-g) zyy=g*cos(kphi(j))**2+d*sin(kphi(j))**2 zyx=zxy zzz=d else ku=kr*sqrt((cos(kphi(i))-cos(kphi(j)))**2 + & (sin(kphi(i))-sin(kphi(j)))**2) a=(1/ku + ii/ku**2 - 1/ku**3) b=(1/ku + ii*3.0/ku**2 - 3.0/ku**3) c=(e-1.0)*ka**2*kdel*kr*cexp(ii*ku)/2.0 zxx=-1.0*c*(a-b*(cos(kphi(i))-cos(kphi(j)))**2/ & (((cos(kphi(i))-cos(kphi(j)))**2 + & (sin(kphi(i))-sin(kphi(j)))**2))) zxy=c*b*(sin(kphi(i))-sin(kphi(j)))* & (cos(kphi(i))-cos(kphi(j)))/ & (((cos(kphi(i))-cos(kphi(j)))**2 + & (sin(kphi(i))-sin(kphi(j)))**2)) zyx=zxy zyy=-1.0*c*(a-b*(sin(kphi(i))-sin(kphi(j)))**2/ & (((cos(kphi(i))-cos(kphi(j)))**2 + & (sin(kphi(i))-sin(kphi(j)))**2))) zzz=-1.0*a*c endif z(3*i-2,3*j-2)=zxx z(3*i-2,3*j-1)=zxy z(3*i-2,3*j)=cmplx(0.0,0.0) z(3*i-1,3*j-2)=zyx z(3*i-1,3*j-1)=zyy z(3*i-1,3*j)=cmplx(0.0,0.0) z(3*i,3*j-2)=cmplx(0.0,0.0) z(3*i,3*j-1)=cmplx(0.0,0.0) z(3*i,3*j)=zzz if(i .eq. 75)then print*, j,ku,zyy,zxy,g,kr,kdel,kl else continue endif 20 continue 10 continue call cgeco(z,dima,orda,piv,reca,wrka) do 26 i=1,3*n pivt(i)=piv(i) do 27 j=1,3*n zt(i,j)=z(i,j) 27 continue 26 continue call cgesl(z,dima,orda,piv,eih,0) call cgesl(zt,dimt,ordt,pivt,eiv,0) open(unit=20,file='circ.out') eaha=0.0 ebha=0.0 echa=0.0 eava=0.0 ebva=0.0 ecva=0.0 print*, 'i,eih(3*i-2),eih(3*i-1),eih(3*i), & eiv(3*i-2),eiv(3*i-1),eiv(3*i)' do 30 i=1,n print*, i,eih(3*i-2),eih(3*i-1),eih(3*i), & eiv(3*i-2),eiv(3*i-1),eiv(3*i) eah=(eih(3*i-2)*cos(kphi(i))+eih(3*i-1)* & sin(kphi(i)))*(e+1.0)/ & (2.0*sin(kphi(i)-phi*degrees)*phase(i)) ebh=(eih(3*i-1)*cos(kphi(i))-eih(3*i-2)*sin(kphi(i)))/ & (cos(kphi(i)-phi*degrees)*phase(i)) ech=eih(3*i) eav=(eiv(3*i-2)*cos(kphi(i))+eiv(3*i-1)* & sin(kphi(i)))*(e+1.0)/ & (2.0*cos(kphi(i)-phi*degrees)* & cos(theta*degrees)*phase(i)) ebv=(eiv(3*i-1)*cos(kphi(i))-eiv(3*i-2)*sin(kphi(i)))/ & (-1.0*sin(kphi(i)-phi*degrees)* & cos(theta*degrees)*phase(i)) ecv=eiv(3*i)*(e+1.0)/(-2.0*sin(theta*degrees)*phase(i)) write(20,*)i,kphi(i),cabs(eah),cabs(ebh),cabs(ech) & ,cabs(eav),cabs(ebv),cabs(ecv) eaha=eaha+cabs(eah)/n ebha=ebha+cabs(ebh)/n echa=echa+cabs(ech)/n eava=eava+cabs(eav)/n ebva=ebva+cabs(ebv)/n ecva=ecva+cabs(ecv)/n c maga=cabs(eiv(i)) c phsa=atan2(real(eiv(i)),aimag(eiv(i))) c magt=cabs(eih(i)) c phst=atan2(real(eih(i)),aimag(eih(i))) 30 continue print*, 'eaha,ebha,echa,eava,ebva,ecva' print*, eaha,ebha,echa,eava,ebva,ecva 100 continue close(20) 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