REAL FUNCTION SCSETRHSSTOP(B,WRK,EPSILON,IPAR,PRECONL,PSCNRM) IMPLICIT NONE * .. Scalar Arguments .. REAL EPSILON * .. * .. Array Arguments .. COMPLEX B(*),WRK(*) INTEGER IPAR(*) * .. * .. Function Arguments .. REAL PSCNRM EXTERNAL PSCNRM * .. * .. Subroutine Arguments .. EXTERNAL PRECONL * .. * .. Local Scalars .. INTEGER LOCLEN,STOPTYPE * .. LOCLEN = IPAR(4) STOPTYPE = IPAR(9) IF ((STOPTYPE.EQ.1) .OR. (STOPTYPE.EQ.4) .OR. + (STOPTYPE.EQ.7)) THEN * ||r||CG iter=',I4,' f.err=',1PE11.3) RETURN END SUBROUTINE PROGRESS(LOCLEN,ITNO,NORMRES,X,RES,TRUERES) C part of PIM C history C 95.08.14 (BTD): modified to include COMMON/NORMERR/ERRSCAL to allow C desired normalization of error. C 96.11.21 (BTD): added IDVOUT2 to COMMON/NORMERR/ to allow passing C of device number for standard output (so that it C doesn't need to be hardwired here). C end history C C Note: when used with STOPTYPE=5, NORMRES=|Ax-b|. C If quantity ERRSCAL is set to |b| elsewhere, then NORMRES/ERRSCAL C is a measure of the fractional error in the solution vector x. IMPLICIT NONE C Arguments: REAL NORMRES INTEGER ITNO,LOCLEN COMPLEX RES(*),TRUERES(*),X(*) C C Common: INTEGER IDVOUT2 REAL ERRSCAL C----------------------------------------------------------------------- COMMON/NORMERR/ERRSCAL,IDVOUT2 C----------------------------------------------------------------------- WRITE(IDVOUT2,FMT=9000)ITNO,NORMRES/ERRSCAL 9000 FORMAT(1X,'iter= ',I4,' frac.err= ',0PF11.7) RETURN END SUBROUTINE CVPROD(N,CX,INCX,CY,INCY) C part of PIM IMPLICIT NONE * * Modified from saxpy level 1 BLAS * element-wise vector multiplication, y<-x*y * Rudnei Dias da Cunha, 16/6/93 * * constant times a vector plus a vector. * uses unrolled loops for increments equal to one. * jack dongarra, linpack, 3/11/78. * * * .. Scalar Arguments .. INTEGER INCX,INCY,N * .. * .. Array Arguments .. COMPLEX CX(*),CY(*) * .. * .. Local Scalars .. INTEGER I,IX,IY,M,MP1 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * 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)*CX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN * * code for both increments equal to 1 * * * clean-up loop * 20 M = MOD(N,4) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M CY(I) = CY(I)*CX(I) 30 CONTINUE IF (N.LT.4) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 CY(I) = CY(I)*CX(I) CY(I+1) = CY(I+1)*CX(I+1) CY(I+2) = CY(I+2)*CX(I+2) CY(I+3) = CY(I+3)*CX(I+3) 50 CONTINUE RETURN END SUBROUTINE PRINTV(N,U) C part of PIM IMPLICIT NONE * .. Scalar Arguments .. INTEGER N * .. * .. Array Arguments .. COMPLEX U(*) * .. * .. Local Scalars .. INTEGER I * .. WRITE (6,FMT=9000) (U(I),I=1,N) RETURN 9000 FORMAT (8 (E14.8,1X)) END SUBROUTINE CINIT(N,ALPHA,CX,INCX) C part of PIM IMPLICIT NONE * Initialises a vector x with a scalar alpha. * Modified from scopy, BLAS Level 1. * Rudnei Dias da Cunha, 14/6/93. * * copies a vector, x, to a vector, y. * uses unrolled loops for increments equal to one. * jack dongarra, linpack, 3/11/78. * * * .. Scalar Arguments .. COMPLEX ALPHA INTEGER INCX,N * .. * .. Array Arguments .. COMPLEX CX(*) * .. * .. Local Scalars .. INTEGER I,IX,M,MP1 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. IF (N.LE.0) RETURN IF (INCX.EQ.1) GO TO 20 * * code for unequal increments or equal increments * not equal to 1 * IX = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 DO 10 I = 1,N CX(IX) = ALPHA IX = IX + INCX 10 CONTINUE RETURN * * code for both increments equal to 1 * * * clean-up loop * 20 M = MOD(N,7) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M CX(I) = ALPHA 30 CONTINUE IF (N.LT.7) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 CX(I) = ALPHA CX(I+1) = ALPHA CX(I+2) = ALPHA CX(I+3) = ALPHA CX(I+4) = ALPHA CX(I+5) = ALPHA CX(I+6) = ALPHA 50 CONTINUE RETURN END COMPLEX FUNCTION CSIGN(X) C part of PIM IMPLICIT NONE * .. Scalar Arguments .. COMPLEX X * .. * .. Intrinsic Functions .. INTRINSIC ABS * .. CSIGN = X/ABS(X) RETURN END SUBROUTINE DECODE(RHO,C,S) C part of PIM IMPLICIT NONE * .. Scalar Arguments .. COMPLEX C,RHO,S * .. * .. Intrinsic Functions .. INTRINSIC ABS,SQRT * .. * .. Parameters .. REAL ONE PARAMETER (ONE=1.0) COMPLEX CZERO PARAMETER (CZERO= (0.0,0.0)) COMPLEX CONE PARAMETER (CONE= (1.0,0.0)) * .. IF (RHO.EQ.CONE) THEN C = CZERO S = CONE ELSE IF (ABS(RHO).LT.ONE) THEN S = 2.0*RHO C = SQRT(CONE-S**2) ELSE C = 2.0/RHO S = SQRT(CONE-C**2) END IF RETURN END SUBROUTINE ENCODE(RHO,C,S) C part of PIM IMPLICIT NONE * .. Scalar Arguments .. COMPLEX C,RHO,S * .. * .. External Functions .. COMPLEX CSIGN EXTERNAL CSIGN * .. * .. Intrinsic Functions .. INTRINSIC ABS * .. * .. Parameters .. COMPLEX CZERO PARAMETER (CZERO= (0.0,0.0)) COMPLEX CONE PARAMETER (CONE= (1.0,0.0)) * .. IF (C.EQ.CZERO) THEN RHO = CONE ELSE IF (ABS(S).LT.ABS(C)) THEN RHO = CSIGN(C)*S/2.0 ELSE RHO = 2.0*CSIGN(S)/C END IF RETURN END SUBROUTINE GIVENS(A,B,C,S) C part of PIM IMPLICIT NONE * .. Scalar Arguments .. COMPLEX A,B,C,S * .. * .. Local Scalars .. COMPLEX TAU * .. * .. Intrinsic Functions .. INTRINSIC ABS,SQRT * .. * .. Parameters .. COMPLEX CZERO PARAMETER (CZERO= (0.0,0.0)) COMPLEX CONE PARAMETER (CONE= (1.0,0.0)) * .. IF (B.EQ.CZERO) THEN C = CONE S = CZERO ELSE IF (ABS(B).GT.ABS(A)) THEN TAU = -A/B S = CONE/SQRT(CONE+TAU**2) C = S*TAU ELSE TAU = -B/A C = CONE/SQRT(CONE+TAU**2) S = C*TAU END IF RETURN END REAL FUNCTION PSCNRM2(LOCLEN,U) C part of PIM package * .. Scalar Arguments .. INTEGER LOCLEN * .. * .. Array Arguments .. COMPLEX U(*) * .. * .. External Functions .. REAL SCNRM2 EXTERNAL SCNRM2 * .. PSCNRM2 = SCNRM2(LOCLEN,U,1) RETURN END SUBROUTINE PCSUM(ISIZE,X) * .. Scalar Arguments .. INTEGER ISIZE * .. * .. Array Arguments .. COMPLEX X(*) * .. RETURN END 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, 4/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 a vector. c jack dongarra, 3/11/78. c complex cx(1),cy(1),ctemp integer i,ix,iy,n,incx,incy ctemp = (0.0e0,0.0e0) cdotc = (0.0e0,0.0e0) 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 a vector. c jack dongarra, 3/11/78. c complex cx(1),cy(1),ctemp integer i,ix,iy,n,incx,incy ctemp = (0.0e0,0.0e0) cdotu = (0.0e0,0.0e0) 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 subroutine crotg(ca,cb,c,s) complex ca,cb,s real c real norm,scale complex alpha if (cabs(ca) .ne. 0.0e0) go to 10 c = 0.0e0 s = (1.0e0,0.0e0) ca = cb go to 20 10 continue scale = cabs(ca) + cabs(cb) norm = scale*sqrt((cabs(ca/cmplx(scale,0.0e0)))**2 + * (cabs(cb/cmplx(scale,0.0e0)))**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, 3/11/78. c modified to correct problem with negative increment, 8/21/90. c complex ca,cx(1) integer i,incx,ix,n c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n cx(ix) = ca*cx(ix) ix = ix + incx 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 c real 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 cswap (n,cx,incx,cy,incy) c c interchanges two vectors. c jack dongarra, 3/11/78. c complex cx(1),cy(1),ctemp integer i,ix,iy,n,incx,incy 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 real function scnrm2( n, cx, incx) logical imag, scale integer i, incx, ix, n, next real cutlo, cuthi, hitest, sum, xmax, absx, zero, one complex cx(1) real real,aimag complex zdumr,zdumi real(zdumr) = zdumr aimag(zdumi) = (0.0e0,-1.0e0)*zdumi 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 modified to correct problem with negative increment, 8/21/90. 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 scnrm2 = zero go to 300 c 10 assign 30 to next sum = zero i = 1 if( incx .lt. 0 )i = (-n+1)*incx + 1 c begin main loop do 220 ix = 1,n 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( 2*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 i = i + incx 220 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(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c modified 9/27/86. c real da,db,c,s,roe,scale,r,z c roe = db if( abs(da) .gt. abs(db) ) roe = da scale = abs(da) + abs(db) if( scale .ne. 0.0e0 ) go to 10 c = 1.0e0 s = 0.0e0 r = 0.0e0 go to 20 10 r = scale*sqrt((da/scale)**2 + (db/scale)**2) r = sign(1.0e0,roe)*r c = da/r s = db/r 20 z = s if( abs(c) .gt. 0.0e0 .and. abs(c) .le. s ) z = 1.0e0/c da = r db = z return end SUBROUTINE DMACHCONS(WHAT,RESULT) * These values are for IEEE-754 arithmetic * .. Parameters .. DOUBLE PRECISION MACHEPS PARAMETER (MACHEPS=2.2204460492503D-16) DOUBLE PRECISION UNDERFLOW PARAMETER (UNDERFLOW=0.22250739D-307) DOUBLE PRECISION OVERFLOW PARAMETER (OVERFLOW=1.7976313D+308) * .. * .. Scalar Arguments .. DOUBLE PRECISION RESULT CHARACTER WHAT * .. IF ((WHAT.EQ.'M') .OR. (WHAT.EQ.'m')) THEN RESULT = MACHEPS ELSE IF ((WHAT.EQ.'U') .OR. (WHAT.EQ.'u')) THEN RESULT = UNDERFLOW ELSE IF ((WHAT.EQ.'O') .OR. (WHAT.EQ.'o')) THEN RESULT = OVERFLOW END IF RETURN END SUBROUTINE PIMDGETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM, + NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT, + ITNO,STATUS,STEPERR,EPSILON,EXITNORM) c PIM -- The Parallel Iterative Methods package c --------------------------------------------- c c Rudnei Dias da Cunha c Centro de Processamento de Dados, c Universidade Federal do Rio Grande do Sul, Brasil c and c Computing Laboratory, University of Kent at Canterbury, U.K. c c Tim Hopkins c Computing Laboratory, University of Kent at Canterbury, U.K. c c ---------------------------------------------------------------------- c c Description of parameter arrays c IPAR (INPUT) : integer c ipar( 1): lda (Leading dimension of a) c 2 : n (Number of rows/columns of a) c 3 : blksz (Size of block of data; used when data is c partitioned using cyclic mode) c 4 : loclen (Number of elements stored locally; c *PARALLEL: Equal to at least m/nprocs or c n/procs depending if row or c column partitioning is used or, c in the case of cyclic partitioning, c it is a multiple of either c m/(blksz*nprocs) or n/(blksz*nprocs). c *SEQUENTIAL: equal to n) c 5 : basisdim (Dimension of orthogonal basis, used in c GMRES) c 6 : nprocs (Number of processors) c 7 : procid (Processor identification) c 8 : precontype (Type of preconditioning; one of c 0 : No preconditioning, c 1 : Left preconditioning, c 2 : Right preconditioning, c 3 : Symmetric preconditioning) c 9 : stoptype (Type of stopping criteria used) c 10 : maxit (Maximum number of iterations allowed) c c IPAR (OUTPUT) : integer c ipar(11): itno (Number of iterations executed) c 12 : status (On exit of iterative method, one of c 0: converged to solution c -1: no convergence has been achieved c -2: "soft"-breakdown, solution may have c been found c -3: "hard"-breakdown, no solution) c -4: conflict in preconditioner and stopping c criterion selected c -5: error in stopping criterion 3, r^{T}z<0) c 13 : steperr (If status is either -2 or -3, it gives c the step number of the respective algorithm c where a breakdown has occurred. Refer to the c User's Guide for further information) c c RPAR/DPAR (INPUT) : real/double precision c dpar( 1): epsilon (The value of epsilon for use in the c stopping criterion) c c RPAR/DPAR (OUTPUT) : real/double precision c dpar( 2): exitnorm (The value of a norm of either the residual c vector or the difference between two c successive solution estimates according to c the value of stoptype) c 3, c 4 : smallest and largest eigenvalues of Q1AQ2 (in the c symmetric case) OR smallest and largest real parts c (in the nonsymmetric case) c 5, c 6 : smallest and largest imaginary parts (only in the c nonsymmetric case) c C .. Parameters .. INTEGER IPARSIZ PARAMETER (IPARSIZ=13) INTEGER DPARSIZ PARAMETER (DPARSIZ=6) C .. C .. Scalar Arguments .. DOUBLE PRECISION EPSILON,EXITNORM INTEGER BASISDIM,BLKSZ,ITNO,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE, + PROCID,STATUS,STEPERR,STOPTYPE C .. C .. Array Arguments .. DOUBLE PRECISION DPAR(DPARSIZ) INTEGER IPAR(IPARSIZ) C .. LDA = IPAR(1) N = IPAR(2) BLKSZ = IPAR(3) LOCLEN = IPAR(4) BASISDIM = IPAR(5) NPROCS = IPAR(6) PROCID = IPAR(7) PRECONTYPE = IPAR(8) STOPTYPE = IPAR(9) MAXIT = IPAR(10) ITNO = IPAR(11) STATUS = IPAR(12) STEPERR = IPAR(13) EPSILON = DPAR(1) EXITNORM = DPAR(2) RETURN END SUBROUTINE PIMDPRTPAR(IPAR,DPAR) c PIM -- The Parallel Iterative Methods package c --------------------------------------------- c c Rudnei Dias da Cunha c Centro de Processamento de Dados, c Universidade Federal do Rio Grande do Sul, Brasil c and c Computing Laboratory, University of Kent at Canterbury, U.K. c c Tim Hopkins c Computing Laboratory, University of Kent at Canterbury, U.K. c c ---------------------------------------------------------------------- c c Description of parameter arrays c IPAR (INPUT) : integer c ipar( 1): lda (Leading dimension of a) c 2 : n (Number of rows/columns of a) c 3 : blksz (Size of block of data; used when data is c partitioned using cyclic mode) c 4 : loclen (Number of elements stored locally; c *PARALLEL: Equal to at least m/nprocs or c n/procs depending if row or c column partitioning is used or, c in the case of cyclic partitioning, c it is a multiple of either c m/(blksz*nprocs) or n/(blksz*nprocs). c *SEQUENTIAL: equal to n) c 5 : basisdim (Dimension of orthogonal basis, used in c GMRES) c 6 : nprocs (Number of processors) c 7 : procid (Processor identification) c 8 : precontype (Type of preconditioning; one of c 0 : No preconditioning, c 1 : Left preconditioning, c 2 : Right preconditioning, c 3 : Symmetric preconditioning) c 9 : stoptype (Type of stopping criteria used) c 10 : maxit (Maximum number of iterations allowed) c c IPAR (OUTPUT) : integer c ipar(11): itno (Number of iterations executed) c 12 : status (On exit of iterative method, one of c 0: converged to solution c -1: no convergence has been achieved c -2: "soft"-breakdown, solution may have c been found c -3: "hard"-breakdown, no solution) c -4: conflict in preconditioner and stopping c criterion selected c -5: error in stopping criterion 3, r^{T}z<0) c 13 : steperr (If status is either -2 or -3, it gives c the step number of the respective algorithm c where a breakdown has occurred. Refer to the c User's Guide for further information) c c RPAR/DPAR (INPUT) : real/double precision c dpar( 1): epsilon (The value of epsilon for use in the c stopping criterion) c c RPAR/DPAR (OUTPUT) : real/double precision c dpar( 2): exitnorm (The value of a norm of either the residual c vector or the difference between two c successive solution estimates according to c the value of stoptype) c 3, c 4 : smallest and largest eigenvalues of Q1AQ2 (in the c symmetric case) OR smallest and largest real parts c (in the nonsymmetric case) c 5, c 6 : smallest and largest imaginary parts (only in the c nonsymmetric case) c C .. Parameters .. INTEGER IPARSIZ PARAMETER (IPARSIZ=13) INTEGER DPARSIZ PARAMETER (DPARSIZ=6) C .. C .. Array Arguments .. DOUBLE PRECISION DPAR(DPARSIZ) INTEGER IPAR(IPARSIZ) C .. C .. Local Scalars .. INTEGER I C .. WRITE (6,FMT=10) (IPAR(I),I=1,IPARSIZ) 10 FORMAT ('lda=',i6,/,'n=',i6,/,'blksz=',i6,/,'loclen=',i4,/, + 'basisdim=',i4,/,'nprocs=',i4,/,'procid=',i4,/, + 'precontype=',i4,/,'stoptype=',i4,/,'maxit=',i4,/,'itno=', + i4,/,'status=',i4,/,'steperr=',i4) WRITE (6,FMT=20) (DPAR(I),I=1,DPARSIZ) 20 FORMAT ('epsilon=',d20.12,/,'exitnorm=',d20.12,/, + 'eigenvalues region:',/,4(d20.12,1x)) RETURN END SUBROUTINE PIMDSETPAR(IPAR,DPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM, + NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT, + EPSILON) c PIM -- The Parallel Iterative Methods package c --------------------------------------------- c c Rudnei Dias da Cunha c Centro de Processamento de Dados, c Universidade Federal do Rio Grande do Sul, Brasil c and c Computing Laboratory, University of Kent at Canterbury, U.K. c c Tim Hopkins c Computing Laboratory, University of Kent at Canterbury, U.K. c c ---------------------------------------------------------------------- c c Description of parameter arrays c IPAR (INPUT) : integer c ipar( 1): lda (Leading dimension of a) c 2 : n (Number of rows/columns of a) c 3 : blksz (Size of block of data; used when data is c partitioned using cyclic mode) c 4 : loclen (Number of elements stored locally; c *PARALLEL: Equal to at least m/nprocs or c n/procs depending if row or c column partitioning is used or, c in the case of cyclic partitioning, c it is a multiple of either c m/(blksz*nprocs) or n/(blksz*nprocs). c *SEQUENTIAL: equal to n) c 5 : basisdim (Dimension of orthogonal basis, used in c GMRES) c 6 : nprocs (Number of processors) c 7 : procid (Processor identification) c 8 : precontype (Type of preconditioning; one of c 0 : No preconditioning, c 1 : Left preconditioning, c 2 : Right preconditioning, c 3 : Symmetric preconditioning) c 9 : stoptype (Type of stopping criteria used) c 10 : maxit (Maximum number of iterations allowed) c c IPAR (OUTPUT) : integer c ipar(11): itno (Number of iterations executed) c 12 : status (On exit of iterative method, one of c 0: converged to solution c -1: no convergence has been achieved c -2: "soft"-breakdown, solution may have c been found c -3: "hard"-breakdown, no solution) c -4: conflict in preconditioner and stopping c criterion selected c -5: error in stopping criterion 3, r^{T}z<0) c 13 : steperr (If status is either -2 or -3, it gives c the step number of the respective algorithm c where a breakdown has occurred. Refer to the c User's Guide for further information) c c RPAR/DPAR (INPUT) : real/double precision c dpar( 1): epsilon (The value of epsilon for use in the c stopping criterion) c c RPAR/DPAR (OUTPUT) : real/double precision c dpar( 2): exitnorm (The value of a norm of either the residual c vector or the difference between two c successive solution estimates according to c the value of stoptype) c 3, c 4 : smallest and largest eigenvalues of Q1AQ2 (in the c symmetric case) OR smallest and largest real parts c (in the nonsymmetric case) c 5, c 6 : smallest and largest imaginary parts (only in the c nonsymmetric case) c C .. Parameters .. DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) INTEGER IPARSIZ PARAMETER (IPARSIZ=13) INTEGER DPARSIZ PARAMETER (DPARSIZ=6) C .. C .. Scalar Arguments .. DOUBLE PRECISION EPSILON INTEGER BASISDIM,BLKSZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE, + PROCID,STOPTYPE C .. C .. Array Arguments .. DOUBLE PRECISION DPAR(DPARSIZ) INTEGER IPAR(IPARSIZ) C .. C .. Local Scalars .. DOUBLE PRECISION EXITNORM INTEGER ITNO,STATUS,STEPERR C .. IPAR(1) = LDA IPAR(2) = N IPAR(3) = BLKSZ IPAR(4) = LOCLEN IPAR(5) = BASISDIM IPAR(6) = NPROCS IPAR(7) = PROCID IPAR(8) = PRECONTYPE IPAR(9) = STOPTYPE IPAR(10) = MAXIT IPAR(11) = -1 IPAR(12) = -1 IPAR(13) = -1 DPAR(1) = EPSILON DPAR(2) = -ONE RETURN END SUBROUTINE PIMSGETPAR(IPAR,SPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM, + NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT, + ITNO,STATUS,STEPERR,EPSILON,EXITNORM) c PIM -- The Parallel Iterative Methods package c --------------------------------------------- c c Rudnei Dias da Cunha c Centro de Processamento de Dados, c Universidade Federal do Rio Grande do Sul, Brasil c and c Computing Laboratory, University of Kent at Canterbury, U.K. c c Tim Hopkins c Computing Laboratory, University of Kent at Canterbury, U.K. c c ---------------------------------------------------------------------- c c Description of parameter arrays c IPAR (INPUT) : integer c ipar( 1): lda (Leading dimension of a) c 2 : n (Number of rows/columns of a) c 3 : blksz (Size of block of data; used when data is c partitioned using cyclic mode) c 4 : loclen (Number of elements stored locally; c *PARALLEL: Equal to at least m/nprocs or c n/procs depending if row or c column partitioning is used or, c in the case of cyclic partitioning, c it is a multiple of either c m/(blksz*nprocs) or n/(blksz*nprocs). c *SEQUENTIAL: equal to n) c 5 : basisdim (Dimension of orthogonal basis, used in c GMRES) c 6 : nprocs (Number of processors) c 7 : procid (Processor identification) c 8 : precontype (Type of preconditioning; one of c 0 : No preconditioning, c 1 : Left preconditioning, c 2 : Right preconditioning, c 3 : Symmetric preconditioning) c 9 : stoptype (Type of stopping criteria used) c 10 : maxit (Maximum number of iterations allowed) c c IPAR (OUTPUT) : integer c ipar(11): itno (Number of iterations executed) c 12 : status (On exit of iterative method, one of c 0: converged to solution c -1: no convergence has been achieved c -2: "soft"-breakdown, solution may have c been found c -3: "hard"-breakdown, no solution) c -4: conflict in preconditioner and stopping c criterion selected c -5: error in stopping criterion 3, r^{T}z<0) c 13 : steperr (If status is either -2 or -3, it gives c the step number of the respective algorithm c where a breakdown has occurred. Refer to the c User's Guide for further information) c c RPAR/DPAR (INPUT) : real/real c spar( 1): epsilon (The value of epsilon for use in the c stopping criterion) c c RPAR/DPAR (OUTPUT) : real/real c spar( 2): exitnorm (The value of a norm of either the residual c vector or the difference between two c successive solution estimates according to c the value of stoptype) c 3, c 4 : smallest and largest eigenvalues of Q1AQ2 (in the c symmetric case) OR smallest and largest real parts c (in the nonsymmetric case) c 5, c 6 : smallest and largest imaginary parts (only in the c nonsymmetric case) c C .. Parameters .. INTEGER IPARSIZ PARAMETER (IPARSIZ=13) INTEGER SPARSIZ PARAMETER (SPARSIZ=6) C .. C .. Scalar Arguments .. REAL EPSILON,EXITNORM INTEGER BASISDIM,BLKSZ,ITNO,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE, + PROCID,STATUS,STEPERR,STOPTYPE C .. C .. Array Arguments .. REAL SPAR(SPARSIZ) INTEGER IPAR(IPARSIZ) C .. LDA = IPAR(1) N = IPAR(2) BLKSZ = IPAR(3) LOCLEN = IPAR(4) BASISDIM = IPAR(5) NPROCS = IPAR(6) PROCID = IPAR(7) PRECONTYPE = IPAR(8) STOPTYPE = IPAR(9) MAXIT = IPAR(10) ITNO = IPAR(11) STATUS = IPAR(12) STEPERR = IPAR(13) EPSILON = SPAR(1) EXITNORM = SPAR(2) RETURN END SUBROUTINE PIMSPRTPAR(IPAR,SPAR) c PIM -- The Parallel Iterative Methods package c --------------------------------------------- c c Rudnei Dias da Cunha c Centro de Processamento de Dados, c Universidade Federal do Rio Grande do Sul, Brasil c and c Computing Laboratory, University of Kent at Canterbury, U.K. c c Tim Hopkins c Computing Laboratory, University of Kent at Canterbury, U.K. c c ---------------------------------------------------------------------- c c Description of parameter arrays c IPAR (INPUT) : integer c ipar( 1): lda (Leading dimension of a) c 2 : n (Number of rows/columns of a) c 3 : blksz (Size of block of data; used when data is c partitioned using cyclic mode) c 4 : loclen (Number of elements stored locally; c *PARALLEL: Equal to at least m/nprocs or c n/procs depending if row or c column partitioning is used or, c in the case of cyclic partitioning, c it is a multiple of either c m/(blksz*nprocs) or n/(blksz*nprocs). c *SEQUENTIAL: equal to n) c 5 : basisdim (Dimension of orthogonal basis, used in c GMRES) c 6 : nprocs (Number of processors) c 7 : procid (Processor identification) c 8 : precontype (Type of preconditioning; one of c 0 : No preconditioning, c 1 : Left preconditioning, c 2 : Right preconditioning, c 3 : Symmetric preconditioning) c 9 : stoptype (Type of stopping criteria used) c 10 : maxit (Maximum number of iterations allowed) c c IPAR (OUTPUT) : integer c ipar(11): itno (Number of iterations executed) c 12 : status (On exit of iterative method, one of c 0: converged to solution c -1: no convergence has been achieved c -2: "soft"-breakdown, solution may have c been found c -3: "hard"-breakdown, no solution) c -4: conflict in preconditioner and stopping c criterion selected c -5: error in stopping criterion 3, r^{T}z<0) c 13 : steperr (If status is either -2 or -3, it gives c the step number of the respective algorithm c where a breakdown has occurred. Refer to the c User's Guide for further information) c c RPAR/DPAR (INPUT) : real/real c spar( 1): epsilon (The value of epsilon for use in the c stopping criterion) c c RPAR/DPAR (OUTPUT) : real/real c spar( 2): exitnorm (The value of a norm of either the residual c vector or the difference between two c successive solution estimates according to c the value of stoptype) c 3, c 4 : smallest and largest eigenvalues of Q1AQ2 (in the c symmetric case) OR smallest and largest real parts c (in the nonsymmetric case) c 5, c 6 : smallest and largest imaginary parts (only in the c nonsymmetric case) c C .. Parameters .. INTEGER IPARSIZ PARAMETER (IPARSIZ=13) INTEGER SPARSIZ PARAMETER (SPARSIZ=6) C .. C .. Array Arguments .. REAL SPAR(SPARSIZ) INTEGER IPAR(IPARSIZ) C .. C .. Local Scalars .. INTEGER I C .. WRITE (6,FMT=10) (IPAR(I),I=1,IPARSIZ) 10 FORMAT ('lda=',i6,/,'n=',i6,/,'blksz=',i6,/,'loclen=',i4,/, + 'basisdim=',i4,/,'nprocs=',i4,/,'procid=',i4,/, + 'precontype=',i4,/,'stoptype=',i4,/,'maxit=',i4,/,'itno=', + i4,/,'status=',i4,/,'steperr=',i4) WRITE (6,FMT=20) (SPAR(I),I=1,SPARSIZ) 20 FORMAT ('epsilon=',e20.12,/,'exitnorm=',e20.12, + 'eigenvalues region:',/,4(e20.12,1x)) RETURN END SUBROUTINE PIMSSETPAR(IPAR,SPAR,LDA,N,BLKSZ,LOCLEN,BASISDIM, + NPROCS,PROCID,PRECONTYPE,STOPTYPE,MAXIT, + EPSILON) c PIM -- The Parallel Iterative Methods package c --------------------------------------------- c c Rudnei Dias da Cunha c Centro de Processamento de Dados, c Universidade Federal do Rio Grande do Sul, Brasil c and c Computing Laboratory, University of Kent at Canterbury, U.K. c c Tim Hopkins c Computing Laboratory, University of Kent at Canterbury, U.K. c c ---------------------------------------------------------------------- c c Description of parameter arrays c IPAR (INPUT) : integer c ipar( 1): lda (Leading dimension of a) c 2 : n (Number of rows/columns of a) c 3 : blksz (Size of block of data; used when data is c partitioned using cyclic mode) c 4 : loclen (Number of elements stored locally; c *PARALLEL: Equal to at least m/nprocs or c n/procs depending if row or c column partitioning is used or, c in the case of cyclic partitioning, c it is a multiple of either c m/(blksz*nprocs) or n/(blksz*nprocs). c *SEQUENTIAL: equal to n) c 5 : basisdim (Dimension of orthogonal basis, used in c GMRES) c 6 : nprocs (Number of processors) c 7 : procid (Processor identification) c 8 : precontype (Type of preconditioning; one of c 0 : No preconditioning, c 1 : Left preconditioning, c 2 : Right preconditioning, c 3 : Symmetric preconditioning) c 9 : stoptype (Type of stopping criteria used) c 10 : maxit (Maximum number of iterations allowed) c c IPAR (OUTPUT) : integer c ipar(11): itno (Number of iterations executed) c 12 : status (On exit of iterative method, one of c 0: converged to solution c -1: no convergence has been achieved c -2: "soft"-breakdown, solution may have c been found c -3: "hard"-breakdown, no solution) c -4: conflict in preconditioner and stopping c criterion selected c -5: error in stopping criterion 3, r^{T}z<0) c 13 : steperr (If status is either -2 or -3, it gives c the step number of the respective algorithm c where a breakdown has occurred. Refer to the c User's Guide for further information) c c RPAR/DPAR (INPUT) : real/real c spar( 1): epsilon (The value of epsilon for use in the c stopping criterion) c c RPAR/DPAR (OUTPUT) : real/real c spar( 2): exitnorm (The value of a norm of either the residual c vector or the difference between two c successive solution estimates according to c the value of stoptype) c 3, c 4 : smallest and largest eigenvalues of Q1AQ2 (in the c symmetric case) OR smallest and largest real parts c (in the nonsymmetric case) c 5, c 6 : smallest and largest imaginary parts (only in the c nonsymmetric case) c C .. Parameters .. REAL ONE PARAMETER (ONE=1.0) INTEGER IPARSIZ PARAMETER (IPARSIZ=13) INTEGER SPARSIZ PARAMETER (SPARSIZ=6) C .. C .. Scalar Arguments .. REAL EPSILON INTEGER BASISDIM,BLKSZ,LDA,LOCLEN,MAXIT,N,NPROCS,PRECONTYPE, + PROCID,STOPTYPE C .. C .. Array Arguments .. REAL SPAR(SPARSIZ) INTEGER IPAR(IPARSIZ) C .. C .. Local Scalars .. REAL EXITNORM INTEGER ITNO,STATUS,STEPERR C .. IPAR(1) = LDA IPAR(2) = N IPAR(3) = BLKSZ IPAR(4) = LOCLEN IPAR(5) = BASISDIM IPAR(6) = NPROCS IPAR(7) = PROCID IPAR(8) = PRECONTYPE IPAR(9) = STOPTYPE IPAR(10) = MAXIT IPAR(11) = -1 IPAR(12) = -1 IPAR(13) = -1 SPAR(1) = EPSILON SPAR(2) = -ONE RETURN END SUBROUTINE SMACHCONS(WHAT,RESULT) * These values are for IEEE-754 arithmetic * .. Parameters .. REAL MACHEPS PARAMETER (MACHEPS=1.19209E-07) REAL UNDERFLOW PARAMETER (UNDERFLOW=0.11754945E-37) REAL OVERFLOW PARAMETER (OVERFLOW=3.40282347E38) * .. * .. Scalar Arguments .. REAL RESULT CHARACTER WHAT * .. IF ((WHAT.EQ.'M') .OR. (WHAT.EQ.'m')) THEN RESULT = MACHEPS ELSE IF ((WHAT.EQ.'U') .OR. (WHAT.EQ.'u')) THEN RESULT = UNDERFLOW ELSE IF ((WHAT.EQ.'O') .OR. (WHAT.EQ.'o')) THEN RESULT = OVERFLOW END IF RETURN END