subroutine caxpy(n,ca,cx,incx,cy,incy) c c constant times a vector plus a vector. c jack dongarra, linpack, 3/11/78. 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)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 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 subroutine ccopy(n,cx,incx,cy,incy) c c copies a vector, x, to a vector, y. c jack dongarra, linpack, 3/11/78. c complex cx(1),cy(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 cy(iy) = 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) = cx(i) 30 continue return end complex function cdotc(n,cx,incx,cy,incy) c c forms the dot product of two vectors, conjugating the first c vector. c jack dongarra, linpack, 3/11/78. 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)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 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 complex function cdotu(n,cx,incx,cy,incy) c c forms the dot product of two vectors. c jack dongarra, linpack, 3/11/78. c complex cx(1),cy(1),ctemp integer i,incx,incy,ix,iy,n c ctemp = (0.0,0.0) cdotu = (0.0,0.0) 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 ctemp = ctemp + cx(ix)*cy(iy) ix = ix + incx iy = iy + incy 10 continue cdotu = ctemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ctemp = ctemp + cx(i)*cy(i) 30 continue cdotu = ctemp return end real function cmach(job) integer job 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 c real eps,tiny,huge,s c eps = 1.0 10 eps = eps/2.0 s = 1.0 + eps if (s .gt. 1.0) go to 10 eps = 2.0*eps cmach =eps if( job .eq. 1) return c s = 1.0 20 tiny = s s = s/16.0 if (s*1.0 .ne. 0.0) go to 20 tiny = (tiny/eps)*100. s = real((1.0,0.0)/cmplx(tiny,0.0)) if (s .ne. 1.0/tiny) tiny = sqrt(tiny) huge = 1.0/tiny if (job .eq. 1) cmach = eps if (job .eq. 2) cmach = tiny if (job .eq. 3) cmach = huge return end subroutine crotg(ca,cb,c,s) complex ca,cb,s real c real norm,scale complex alpha if (cabs(ca) .ne. 0.) go to 10 c = 0. s = (1.,0.) ca = cb go to 20 10 continue scale = cabs(ca) + cabs(cb) norm = scale * sqrt((cabs(ca/scale))**2 + (cabs(cb/scale))**2) alpha = ca /cabs(ca) c = cabs(ca) / norm s = alpha * conjg(cb) / norm ca = alpha * norm 20 continue return end subroutine cscal(n,ca,cx,incx) c c scales a vector by a constant. c jack dongarra, linpack, 3/11/78. c complex ca,cx(1) integer i,incx,n,nincx c if(n.le.0)return if(incx.eq.1)go to 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 subroutine csrot (n,cx,incx,cy,incy,c,s) c c applies a plane rotation, where the cos and sin (c and s) are real c and the vectors cx and cy are complex. c jack dongarra, linpack, 3/11/78. c complex cx(1),cy(1),ctemp real 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 ctemp = c*cx(ix) + s*cy(iy) cy(iy) = c*cy(iy) - s*cx(ix) cx(ix) = ctemp 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 ctemp = c*cx(i) + s*cy(i) cy(i) = c*cy(i) - s*cx(i) cx(i) = ctemp 30 continue return end subroutine csscal(n,sa,cx,incx) c c scales a complex vector by a real constant. c jack dongarra, linpack, 3/11/78. c complex cx(1) real sa integer i,incx,n,nincx c if(n.le.0)return if(incx.eq.1)go to 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 subroutine cswap (n,cx,incx,cy,incy) c c interchanges two vectors. c jack dongarra, linpack, 3/11/78. c complex cx(1),cy(1),ctemp 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 ctemp = cx(ix) cx(ix) = cy(iy) cy(iy) = ctemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 20 do 30 i = 1,n ctemp = cx(i) cx(i) = cy(i) cy(i) = ctemp 30 continue return end integer function icamax(n,cx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. 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 = 0 if( n .lt. 1 ) return icamax = 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 = 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 real function scasum(n,cx,incx) c c takes the sum of the absolute values of a complex vector and c returns a single precision result. c jack dongarra, linpack, 3/11/78. 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)go to 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 real function scnrm2( n, cx, incx) logical imag, scale integer next real cutlo, cuthi, hitest, sum, xmax, absx, zero, one complex cx(1) data zero, one /0.0e0, 1.0e0/ c c unitary norm of the complex n-vector stored in cx() 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 / 4.441e-16, 1.304e19 / c if(n .gt. 0) go to 10 scnrm2 = 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 = abs(real(cx(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 = abs(aimag(cx(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 scnrm2 = sqrt(sum) if(scale) scnrm2 = scnrm2 * xmax 300 continue return end subroutine srotg(sa,sb,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c real sa,sb,c,s,roe,scale,r,z c roe = sb if( abs(sa) .gt. abs(sb) ) roe = sa scale = abs(sa) + abs(sb) if( scale .ne. 0.0 ) go to 10 c = 1.0 s = 0.0 r = 0.0 go to 20 10 r = scale*sqrt((sa/scale)**2 + (sb/scale)**2) r = sign(1.0,roe)*r c = sa/r s = sb/r 20 z = 1.0 if( abs(sa) .gt. abs(sb) ) z = s if( abs(sb) .ge. abs(sa) .and. c .ne. 0.0 ) z = 1.0/c sa = r sb = z return end