c c takes the sum of the absolute values. c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision function dzasum(n,zx,incx) double complex zx(*) double precision stemp,abs integer i,incx,ix,n c dzasum = 0.0d0 stemp = 0.0d0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n stemp = stemp + abs(zx(ix)) ix = ix + incx 10 continue dzasum = stemp return c c code for increment equal to 1 c 20 do 30 i = 1,n stemp = stemp + abs(zx(i)) 30 continue dzasum = stemp return end * * DZNRM2 returns the euclidean norm of a vector via the function * name, so that * * DZNRM2 := sqrt( conjg( x' )*x ) * * * * -- This version written on 25-October-1982. * Modified on 14-October-1993 to inline the call to ZLASSQ. * Sven Hammarling, Nag Ltd. * * DOUBLE PRECISION FUNCTION DZNRM2( N, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N * .. Array Arguments .. COMPLEX*16 X( * ) * .. * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. INTEGER IX DOUBLE PRECISION NORM, SCALE, SSQ, TEMP * .. Intrinsic Functions .. INTRINSIC ABS, DIMAG, DBLE, SQRT * .. * .. Executable Statements .. IF( N.LT.1 .OR. INCX.LT.1 )THEN NORM = ZERO ELSE SCALE = ZERO SSQ = ONE * The following loop is equivalent to this call to the LAPACK * auxiliary routine: * CALL ZLASSQ( N, X, INCX, SCALE, SSQ ) * DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX IF( DBLE( X( IX ) ).NE.ZERO )THEN TEMP = ABS( DBLE( X( IX ) ) ) IF( SCALE.LT.TEMP )THEN SSQ = ONE + SSQ*( SCALE/TEMP )**2 SCALE = TEMP ELSE SSQ = SSQ + ( TEMP/SCALE )**2 END IF END IF IF( DIMAG( X( IX ) ).NE.ZERO )THEN TEMP = ABS( DIMAG( X( IX ) ) ) IF( SCALE.LT.TEMP )THEN SSQ = ONE + SSQ*( SCALE/TEMP )**2 SCALE = TEMP ELSE SSQ = SSQ + ( TEMP/SCALE )**2 END IF END IF 10 CONTINUE NORM = SCALE * SQRT( SSQ ) END IF * DZNRM2 = NORM RETURN END c c finds the index of element having max. absolute value. c jack dongarra, 1/15/85. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c integer function izamax(n,zx,incx) double complex zx(*) double precision smax integer i,incx,ix,n c izamax = 0 if( n.lt.1 .or. incx.le.0 )return izamax = 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(zx(1)) ix = ix + incx do 10 i = 2,n if(abs(zx(ix)).le.smax) go to 5 izamax = i smax = abs(zx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 smax = abs(zx(1)) do 30 i = 2,n if(abs(zx(i)).le.smax) go to 30 izamax = i smax = abs(zx(i)) 30 continue return end c c constant times a vector plus a vector. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c subroutine zaxpy(n,za,zx,incx,zy,incy) double complex zx(*),zy(*),za integer i,incx,incy,ix,iy,n double precision abs if(n.le.0)return if (abs(za) .eq. 0.0d0) 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 zy(iy) = zy(iy) + za*zx(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 zy(i) = zy(i) + za*zx(i) 30 continue return end c c copies a vector, x, to a vector, y. c jack dongarra, linpack, 4/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c subroutine zcopy(n,zx,incx,zy,incy) double complex zx(*),zy(*) 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 zy(iy) = zx(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 zy(i) = zx(i) 30 continue return end c c forms the dot product of a vector. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex function zdotc(n,zx,incx,zy,incy) double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n ztemp = (0.0d0,0.0d0) zdotc = (0.0d0,0.0d0) 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 ztemp = ztemp + conjg(zx(ix))*zy(iy) ix = ix + incx iy = iy + incy 10 continue zdotc = ztemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = ztemp + conjg(zx(i))*zy(i) 30 continue zdotc = ztemp return end c c forms the dot product of two vectors. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex function zdotu(n,zx,incx,zy,incy) double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n ztemp = (0.0d0,0.0d0) zdotu = (0.0d0,0.0d0) 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 ztemp = ztemp + zx(ix)*zy(iy) ix = ix + incx iy = iy + incy 10 continue zdotu = ztemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = ztemp + zx(i)*zy(i) 30 continue zdotu = ztemp return end c c scales a vector by a constant. c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c subroutine zdscal(n,da,zx,incx) double complex zx(*) double precision da integer i,incx,ix,n c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n zx(ix) = dcmplx(da,0.0d0)*zx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 do 30 i = 1,n zx(i) = dcmplx(da,0.0d0)*zx(i) 30 continue return end subroutine zrotg(ca,cb,c,s) double complex ca,cb,s double precision c double precision norm,scale double complex alpha if (abs(ca) .ne. 0.0d0) go to 10 c = 0.0d0 s = (1.0d0,0.0d0) ca = cb go to 20 10 continue scale = abs(ca) + abs(cb) norm = scale*dsqrt((abs(ca/cmplx(scale,0.0d0)))**2 + * (abs(cb/cmplx(scale,0.0d0)))**2) alpha = ca /abs(ca) c = abs(ca) / norm s = alpha * conjg(cb) / norm ca = alpha * norm 20 continue return end c c scales a vector by a constant. c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c subroutine zscal(n,za,zx,incx) double complex za,zx(*) integer i,incx,ix,n c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n zx(ix) = za*zx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 do 30 i = 1,n zx(i) = za*zx(i) 30 continue return end c c interchanges two vectors. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c subroutine zswap (n,zx,incx,zy,incy) double complex zx(*),zy(*),ztemp 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 ztemp = zx(ix) zx(ix) = zy(iy) zy(iy) = ztemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 20 do 30 i = 1,n ztemp = zx(i) zx(i) = zy(i) zy(i) = ztemp 30 continue return end