double precision function dcabs1(z) double complex z,zz double precision t(2) equivalence (zz,t(1)) zz = z dcabs1 = dabs(t(1)) + dabs(t(2)) return end subroutine drotg(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c double precision da,db,c,s,roe,scale,r,z c roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .ne. 0.0d0 ) go to 10 c = 1.0d0 s = 0.0d0 r = 0.0d0 go to 20 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r 20 z = 1.0d0 if( dabs(da) .gt. dabs(db) ) z = s if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c da = r db = z return end double precision function dzasum(n,zx,incx) c c takes the sum of the absolute values. c jack dongarra, 3/11/78. c double complex zx(1) double precision stemp,dcabs1 c dzasum = 0.0d0 stemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increments not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n stemp = stemp + dcabs1(zx(ix)) ix = ix + incx 10 continue dzasum = stemp return c c code for increments equal to 1 c 20 do 30 i = 1,n stemp = stemp + dcabs1(zx(i)) 30 continue dzasum = stemp return end double precision function dznrm2( n, zx, incx) logical imag, scale integer next double precision cutlo, cuthi, hitest, sum, xmax, absx, zero, one double complex zx(1) double precision dreal,dimag double complex zdumr,zdumi dreal(zdumr) = zdumr dimag(zdumi) = (0.0d0,-1.0d0)*zdumi data zero, one /0.0d0, 1.0d0/ c c unitary norm of the complex n-vector stored in zx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson , 1978 jan 08 c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of sqrt(u/eps) over all known machines. c cuthi = minimum of sqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 8.232d-11, 1.304d19 / c if(n .gt. 0) go to 10 dznrm2 = zero go to 300 c 10 assign 30 to next sum = zero nn = n * incx c begin main loop do 210 i=1,nn,incx absx = dabs(dreal(zx(i))) imag = .false. go to next,(30, 50, 70, 90, 110) 30 if( absx .gt. cutlo) go to 85 assign 50 to next scale = .false. c c phase 1. sum is zero c 50 if( absx .eq. zero) go to 200 if( absx .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 assign 110 to next sum = (sum / absx) / absx 105 scale = .true. xmax = absx go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( absx .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( absx .le. xmax ) go to 115 sum = one + sum * (xmax / absx)**2 xmax = absx go to 200 c 115 sum = sum + (absx/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c 85 assign 90 to next scale = .false. c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c 90 if(absx .ge. hitest) go to 100 sum = sum + absx**2 200 continue c control selection of real and imaginary parts. c if(imag) go to 210 absx = dabs(dimag(zx(i))) imag = .true. go to next,( 50, 70, 90, 110 ) c 210 continue c c end of main loop. c compute square root and adjust for scaling. c dznrm2 = dsqrt(sum) if(scale) dznrm2 = dznrm2 * xmax 300 continue return end subroutine zdrot (n,zx,incx,zy,incy,c,s) c c applies a plane rotation, where the cos and sin (c and s) are c double precision and the vectors zx and zy are double complex. c jack dongarra, linpack, 3/11/78. c double complex zx(1),zy(1),ztemp double precision c,s 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 = c*zx(ix) + s*zy(iy) zy(iy) = c*zy(iy) - s*zx(ix) zx(ix) = ztemp 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 ztemp = c*zx(i) + s*zy(i) zy(i) = c*zy(i) - s*zx(i) zx(i) = ztemp 30 continue return end integer function izamax(n,zx,incx) c c finds the index of element having max. absolute value. c jack dongarra, 3/11/78. c double complex zx(1) double precision smax double complex zdum double precision dcabs1 c izamax = 1 if(n.le.1)return if(incx.eq.1)go to 20 c c code for increments not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 smax = dcabs1(zx(1)) ix = ix + incx do 10 i = 2,n if(dcabs1(zx(ix)).le.smax) go to 5 izamax = ix smax = dcabs1(zx(ix)) 5 ix = ix + incx 10 continue return c c code for increments equal to 1 c 20 smax = dcabs1(zx(1)) do 30 i = 2,n if(dcabs1(zx(i)).le.smax) go to 30 izamax = i smax = dcabs1(zx(i)) 30 continue return end subroutine zaxpy(n,za,zx,incx,zy,incy) c c constant times a vector plus a vector. c jack dongarra, 3/11/78. c double complex zx(1),zy(1),za double precision dcabs1 if(n.le.0)return if (dcabs1(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 subroutine zcopy(n,zx,incx,zy,incy) c c copies a vector, x, to a vector, y. c jack dongarra, linpack, 4/11/78. c double complex zx(1),zy(1) 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 double complex function zdotc(n,zx,incx,zy,incy) c c forms the dot product of a vector. c jack dongarra, 3/11/78. c double complex zx(1),zy(1),ztemp 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 + dconjg(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 + dconjg(zx(i))*zy(i) 30 continue zdotc = ztemp return end double complex function zdotu(n,zx,incx,zy,incy) c c forms the dot product of a vector. c jack dongarra, 3/11/78. c double complex zx(1),zy(1),ztemp 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 subroutine zdscal(n,da,zx,incx) c c scales a vector by a constant. c jack dongarra, 3/11/78. c double complex zx(1) double precision da if(n.le.0)return if(incx.eq.1)go to 20 c c code for increments not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n zx(ix) = dcmplx(da,0.0d0)*zx(ix) ix = ix + incx 10 continue return c c code for increments equal to 1 c 20 do 30 i = 1,n zx(i) = dcmplx(da,0.0d0)*zx(i) 30 continue return end double precision function zmach(job) integer job c c double complex floating point arithmetic constants. c for the linpack test drivers only. c not used by actual linpack subroutines. c c smach computes machine parameters of floating point c arithmetic for use in testing only. not required by c linpack proper. c c if trouble with automatic computation of these quantities, c they can be set by direct assignment statements. c assume the computer has c c b = base of arithmetic c t = number of base b digits c l = smallest possible exponent c u = largest possible exponent c c then c c eps = b**(1-t) c tiny = 100.0*b**(-l+t) c huge = 0.01*b**(u-t) c c dmach same as smach except t, l, u apply to c double precision. c c cmach same as smach except if complex division c is done by c c 1/(x+i*y) = (x-i*y)/(x**2+y**2) c c then c c tiny = sqrt(tiny) c huge = sqrt(huge) c c c job is 1, 2 or 3 for epsilon, tiny and huge, respectively. c double precision eps,tiny,huge,s c eps = 1.0d0 10 eps = eps/2.0d0 s = 1.0d0 + eps if (s .gt. 1.0d0) go to 10 eps = 2.0d0*eps c s = 1.0d0 20 tiny = s s = s/16.0d0 if (s*1.0d0 .ne. 0.0d0) go to 20 tiny = tiny/eps s = (1.0d0,0.0d0)/dcmplx(tiny,0.0d0) if (s .ne. 1.0d0/tiny) tiny = dsqrt(tiny) huge = 1.0d0/tiny c if (job .eq. 1) zmach = eps if (job .eq. 2) zmach = tiny if (job .eq. 3) zmach = huge return end subroutine zrotg(ca,cb,c,s) double complex ca,cb,s double precision c double precision norm,scale double complex alpha if (cdabs(ca) .ne. 0.0d0) go to 10 c = 0.0d0 s = (1.0d0,0.0d0) ca = cb go to 20 10 continue scale = cdabs(ca) + cdabs(cb) norm = scale*dsqrt((cdabs(ca/dcmplx(scale,0.0d0)))**2 + * (cdabs(cb/dcmplx(scale,0.0d0)))**2) alpha = ca /cdabs(ca) c = cdabs(ca) / norm s = alpha * dconjg(cb) / norm ca = alpha * norm 20 continue return end subroutine zscal(n,za,zx,incx) c c scales a vector by a constant. c jack dongarra, 3/11/78. c double complex za,zx(1) if(n.le.0)return if(incx.eq.1)go to 20 c c code for increments not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n zx(ix) = za*zx(ix) ix = ix + incx 10 continue return c c code for increments equal to 1 c 20 do 30 i = 1,n zx(i) = za*zx(i) 30 continue return end subroutine zswap (n,zx,incx,zy,incy) c c interchanges two vectors. c jack dongarra, 3/11/78. c double complex zx(1),zy(1),ztemp 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