w3profsmd.F90 File Reference

Go to the source code of this file.

Modules
module	w3profsmd

Functions/Subroutines
subroutine	w3profsmd::w3xypug (ISP, FACX, FACY, DTG, VQ, VGX, VGY, LCALC)
subroutine	w3profsmd::w3cflug (ISEA, NKCFL, FACX, FACY, DT, MAPFS, CFLXYMAX, VGX, VGY)
subroutine	w3profsmd::w3xypfsn2 (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3profsmd::w3xypfspsi2 (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3profsmd::w3xypfsnimp (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3profsmd::w3xypfsfct2 (ISP, C, LCALC, RD10, RD20, DT, AC)
subroutine	w3profsmd::setdepth
subroutine	bcgstab (n, rhs, sol, ipar, fpar, w)
subroutine	implu (np, umm, beta, ypiv, u, permut, full)
subroutine	uppdir (n, p, np, lbp, indp, y, u, usav, flops)
subroutine	givens (x, y, c, s)
logical function	stopbis (n, ipar, mvpi, fpar, r, delx, sx)
subroutine	tidycg (n, ipar, fpar, sol, delx)
logical function	brkdn (alpha, ipar)
subroutine	bisinit (ipar, fpar, wksize, dsc, lp, rp, wk)
subroutine	mgsro (full, lda, n, m, ind, ops, vec, hh, ierr)
subroutine	amux (n, x, y, a, ja, ia)
subroutine	amuxms (n, x, y, a, ja)
subroutine	atmux (n, x, y, a, ja, ia)
subroutine	atmuxr (m, n, x, y, a, ja, ia)
subroutine	amuxe (n, x, y, na, ncol, a, ja)
subroutine	amuxd (n, x, y, diag, ndiag, idiag, ioff)
subroutine	amuxj (n, x, y, jdiag, a, ja, ia)
subroutine	vbrmv (nr, nc, ia, ja, ka, a, kvstr, kvstc, x, b)
subroutine	lsol (n, x, y, al, jal, ial)
subroutine	ldsol (n, x, y, al, jal)
subroutine	lsolc (n, x, y, al, jal, ial)
subroutine	ldsolc (n, x, y, al, jal)
subroutine	ldsoll (n, x, y, al, jal, nlev, lev, ilev)
subroutine	usol (n, x, y, au, jau, iau)
subroutine	udsol (n, x, y, au, jau)
subroutine	usolc (n, x, y, au, jau, iau)
subroutine	udsolc (n, x, y, au, jau)
subroutine	lusol (n, y, x, alu, jlu, ju)
subroutine	lutsol (n, y, x, alu, jlu, ju)
subroutine	qsplit (a, ind, n, ncut)
subroutine	runrc (n, rhs, sol, ipar, fpar, wk, guess, a, ja, ia, au, jau, ju, solver)
subroutine	ilut (n, a, ja, ia, lfil, droptol, alu, jlu, ju, iwk, w, jw, ierr)
subroutine	ilu0 (n, a, ja, ia, alu, jlu, ju, iw, ierr)
subroutine	pgmres (n, im, rhs, sol, eps, maxits, aspar, nnz, ia, ja, alu, jlu, ju, vv, ierr)
double precision function	dnrm2 (N, X)
subroutine	dlassq (N, X, SCALE, SUMSQ)
double precision function	ddot (n, dx, dy)
subroutine	daxpy (n, da, dx, incx, dy, incy)

Function/Subroutine Documentation

amux()

subroutine amux	(	integer	n,
		real*8, dimension(*)	x,
		real*8, dimension(*)	y,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(*)	ia
	)

Definition at line 2820 of file w3profsmd.F90.

  real*8  x(*), y(*), a(*)
  integer n, ja(*), ia(*)
  !-----------------------------------------------------------------------
  !         A times a vector
  !-----------------------------------------------------------------------
  ! multiplies a matrix by a vector using the dot product form
  ! Matrix A is stored in compressed sparse row storage.
  !
  ! on entry:
  !----------
  ! n     = row dimension of A
  ! x     = real array of length equal to the column dimension of
  !         the A matrix.
  ! a, ja,
  !    ia = input matrix in compressed sparse row format.
  !
  ! on return:
  !-----------
  ! y     = real array of length n, containing the product y=Ax
  !
  !-----------------------------------------------------------------------
  ! local variables
  !
  real*8 t
  integer i, k
  !-----------------------------------------------------------------------
  do i = 1,n
    !
    !     compute the inner product of row i with vector x
    !
    t = 0.0d0
    do k=ia(i), ia(i+1)-1
      t = t + a(k)*x(ja(k))
    end do
    !
    !     store result in y(i)
    !
    y(i) = t
  end do
  !
  return
  !---------end-of-amux---------------------------------------------------
  !-----------------------------------------------------------------------

Referenced by pgmres(), and runrc().

amuxd()

subroutine amuxd	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(ndiag,idiag)	diag,
		integer	ndiag,
		integer	idiag,
		integer, dimension(idiag)	ioff
	)

Definition at line 3055 of file w3profsmd.F90.

  integer n, ndiag, idiag, ioff(idiag)
  real*8 x(n), y(n), diag(ndiag,idiag)
  !-----------------------------------------------------------------------
  !        A times a vector in Diagonal storage format (DIA)
  !-----------------------------------------------------------------------
  ! multiplies a matrix by a vector when the original matrix is stored
  ! in the diagonal storage format.
  !-----------------------------------------------------------------------
  !
  ! on entry:
  !----------
  ! n     = row dimension of A
  ! x     = real array of length equal to the column dimension of
  !         the A matrix.
  ! ndiag  = integer. The first dimension of array adiag as declared in
  !         the calling program.
  ! idiag  = integer. The number of diagonals in the matrix.
  ! diag   = real array containing the diagonals stored of A.
  ! idiag  = number of diagonals in matrix.
  ! diag   = real array of size (ndiag x idiag) containing the diagonals
  !
  ! ioff   = integer array of length idiag, containing the offsets of the
  !          diagonals of the matrix:
  !          diag(i,k) contains the element a(i,i+ioff(k)) of the matrix.
  !
  ! on return:
  !-----------
  ! y     = real array of length n, containing the product y=A*x
  !
  !-----------------------------------------------------------------------
  !       local variables
  !
  integer j, k, io, i1, i2
  !-----------------------------------------------------------------------
  do j=1, n
    y(j) = 0.0d0
  end do
  do j=1, idiag
    io = ioff(j)
    i1 = max0(1,1-io)
    i2 = min0(n,n-io)
    do k=i1, i2
      y(k) = y(k)+diag(k,j)*x(k+io)
    end do
  end do
  !
  return
  !----------end-of-amuxd-------------------------------------------------
  !-----------------------------------------------------------------------

amuxe()

subroutine amuxe	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		integer	na,
		integer	ncol,
		real*8, dimension(na,*)	a,
		integer, dimension(na,*)	ja
	)

Definition at line 3006 of file w3profsmd.F90.

  implicit none
 
  integer     :: n, na, ncol, ja(na,*)
  real*8      :: x(n), y(n), a(na,*)
 
  !-----------------------------------------------------------------------
  !        A times a vector in Ellpack Itpack format (ELL)
  !-----------------------------------------------------------------------
  ! multiplies a matrix by a vector when the original matrix is stored
  ! in the ellpack-itpack sparse format.
  !-----------------------------------------------------------------------
  !
  ! on entry:
  !----------
  ! n     = row dimension of A
  ! x     = real array of length equal to the column dimension of
  !         the A matrix.
  ! na    = integer. The first dimension of arrays a and ja
  !         as declared by the calling program.
  ! ncol  = integer. The number of active columns in array a.
  !         (i.e., the number of generalized diagonals in matrix.)
  ! a, ja = the real and integer arrays of the itpack format
  !         (a(i,k),k=1,ncol contains the elements of row i in matrix
  !          ja(i,k),k=1,ncol contains their column numbers)
  !
  ! on return:
  !-----------
  ! y     = real array of length n, containing the product y=y=A*x
  !
  !-----------------------------------------------------------------------
  ! local variables
  !
  integer i, j
  !-----------------------------------------------------------------------
  do i=1, n
    y(i) = 0.0
  end do
  do j=1,ncol
    do i = 1,n
      y(i) = y(i)+a(i,j)*x(ja(i,j))
    end do
  end do
  !
  return
  !--------end-of-amuxe---------------------------------------------------
  !-----------------------------------------------------------------------

amuxj()

subroutine amuxj	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		integer	jdiag,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(*)	ia
	)

Definition at line 3107 of file w3profsmd.F90.

  integer n, jdiag, ja(*), ia(*)
  real*8 x(n), y(n), a(*)
  !-----------------------------------------------------------------------
  !        A times a vector in Jagged-Diagonal storage format (JAD)
  !-----------------------------------------------------------------------
  ! multiplies a matrix by a vector when the original matrix is stored
  ! in the jagged diagonal storage format.
  !-----------------------------------------------------------------------
  !
  ! on entry:
  !----------
  ! n      = row dimension of A
  ! x      = real array of length equal to the column dimension of
  !         the A matrix.
  ! jdiag  = integer. The number of jadded-diagonals in the data-structure.
  ! a      = real array containing the jadded diagonals of A stored
  !          in succession (in decreasing lengths)
  ! j      = integer array containing the colum indices of the
  !          corresponding elements in a.
  ! ia     = integer array containing the lengths of the  jagged diagonals
  !
  ! on return:
  !-----------
  ! y      = real array of length n, containing the product y=A*x
  !
  ! Note:
  !-------
  ! Permutation related to the JAD format is not performed.
  ! this can be done by:
  !     call permvec (n,y,y,iperm)
  ! after the call to amuxj, where iperm is the permutation produced
  ! by csrjad.
  !-----------------------------------------------------------------------
  ! local variables
  !
  integer i, ii, k1, ilen, j
  !-----------------------------------------------------------------------
  do i=1, n
    y(i) = 0.0d0
  end do
  do  ii=1, jdiag
    k1 = ia(ii)-1
    ilen = ia(ii+1)-k1-1
    do  j=1,ilen
      y(j)= y(j)+a(k1+j)*x(ja(k1+j))
    end do
  end do
  !
  return
  !----------end-of-amuxj-------------------------------------------------
  !-----------------------------------------------------------------------

amuxms()

subroutine amuxms	(	integer	n,
		real*8, dimension(*)	x,
		real*8, dimension(*)	y,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja
	)

Definition at line 2866 of file w3profsmd.F90.

  real*8  x(*), y(*), a(*)
  integer n, ja(*)
  !-----------------------------------------------------------------------
  !         A times a vector in MSR format
  !-----------------------------------------------------------------------
  ! multiplies a matrix by a vector using the dot product form
  ! Matrix A is stored in Modified Sparse Row storage.
  !
  ! on entry:
  !----------
  ! n     = row dimension of A
  ! x     = real array of length equal to the column dimension of
  !         the A matrix.
  ! a, ja,= input matrix in modified compressed sparse row format.
  !
  ! on return:
  !-----------
  ! y     = real array of length n, containing the product y=Ax
  !
  !-----------------------------------------------------------------------
  ! local variables
  !
  integer i, k
  !-----------------------------------------------------------------------
  do i=1, n
    y(i) = a(i)*x(i)
  end do
  do i = 1,n
    !
    !     compute the inner product of row i with vector x
    !
    do k=ja(i), ja(i+1)-1
      y(i) = y(i) + a(k) *x(ja(k))
    end do
  end do
  !
  return
  !---------end-of-amuxm--------------------------------------------------
  !-----------------------------------------------------------------------

atmux()

subroutine atmux	(	integer	n,
		real*8, dimension(*)	x,
		real*8, dimension(*)	y,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(*)	ia
	)

Definition at line 2908 of file w3profsmd.F90.

  real*8 x(*), y(*), a(*)
  integer n, ia(*), ja(*)
  !-----------------------------------------------------------------------
  !         transp( A ) times a vector
  !-----------------------------------------------------------------------
  ! multiplies the transpose of a matrix by a vector when the original
  ! matrix is stored in compressed sparse row storage. Can also be
  ! viewed as the product of a matrix by a vector when the original
  ! matrix is stored in the compressed sparse column format.
  !-----------------------------------------------------------------------
  !
  ! on entry:
  !----------
  ! n     = row dimension of A
  ! x     = real array of length equal to the column dimension of
  !         the A matrix.
  ! a, ja,
  !    ia = input matrix in compressed sparse row format.
  !
  ! on return:
  !-----------
  ! y     = real array of length n, containing the product y=transp(A)*x
  !
  !-----------------------------------------------------------------------
  !     local variables
  !
  integer i, k
  !-----------------------------------------------------------------------
  !
  !     zero out output vector
  !
  do i=1,n
    y(i) = 0.0
  end do
  !
  ! loop over the rows
  !
  do i = 1,n
    do  k=ia(i), ia(i+1)-1
      y(ja(k)) = y(ja(k)) + x(i)*a(k)
    end do
  end do
  !
  return
  !-------------end-of-atmux----------------------------------------------
  !-----------------------------------------------------------------------

Referenced by runrc().

atmuxr()

subroutine atmuxr	(	integer	m,
		integer	n,
		real*8, dimension(*)	x,
		real*8, dimension(*)	y,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(*)	ia
	)

Definition at line 2957 of file w3profsmd.F90.

  real*8 x(*), y(*), a(*)
  integer m, n, ia(*), ja(*)
  !-----------------------------------------------------------------------
  !         transp( A ) times a vector, A can be rectangular
  !-----------------------------------------------------------------------
  ! See also atmux.  The essential difference is how the solution vector
  ! is initially zeroed.  If using this to multiply rectangular CSC
  ! matrices by a vector, m number of rows, n is number of columns.
  !-----------------------------------------------------------------------
  !
  ! on entry:
  !----------
  ! m     = column dimension of A
  ! n     = row dimension of A
  ! x     = real array of length equal to the column dimension of
  !         the A matrix.
  ! a, ja,
  !    ia = input matrix in compressed sparse row format.
  !
  ! on return:
  !-----------
  ! y     = real array of length n, containing the product y=transp(A)*x
  !
  !-----------------------------------------------------------------------
  !     local variables
  !
  integer i, k
  !-----------------------------------------------------------------------
  !
  !     zero out output vector
  !
  do i=1,m
    y(i) = 0.0
  end do
  !
  ! loop over the rows
  !
  do i = 1,n
    do k=ia(i), ia(i+1)-1
      y(ja(k)) = y(ja(k)) + x(i)*a(k)
    end do
  end do
  !
  return
  !-------------end-of-atmuxr---------------------------------------------
  !-----------------------------------------------------------------------

bcgstab()

subroutine bcgstab	(	integer	n,
		real*8, dimension(n)	rhs,
		real*8, dimension(n)	sol,
		integer, dimension(16)	ipar,
		real*8, dimension(16)	fpar,
		real*8, dimension(n,8)	w
	)

Definition at line 2044 of file w3profsmd.F90.

  implicit none
  integer n, ipar(16)
  real*8 rhs(n), sol(n), fpar(16), w(n,8)
  !-----------------------------------------------------------------------
  !     BCGSTAB --- Bi Conjugate Gradient stabilized (BCGSTAB)
  !     This is an improved BCG routine. (1) no matrix transpose is
  !     involved. (2) the convergence is smoother.
  !
  !
  !     Algorithm:
  !     Initialization - r = b - A x, r0 = r, p = r, rho = (r0, r),
  !     Iterate -
  !     (1) v = A p
  !     (2) alpha = rho / (r0, v)
  !     (3) s = r - alpha v
  !     (4) t = A s
  !     (5) omega = (t, s) / (t, t)
  !     (6) x = x + alpha * p + omega * s
  !     (7) r = s - omega * t
  !     convergence test goes here
  !     (8) beta = rho, rho = (r0, r), beta = rho * alpha / (beta * omega)
  !         p = r + beta * (p - omega * v)
  !
  !     in this routine, before successful return, the fpar's are
  !     fpar(3) == initial (preconditionied-)residual norm
  !     fpar(4) == target (preconditionied-)residual norm
  !     fpar(5) == current (preconditionied-)residual norm
  !     fpar(6) == current residual norm or error
  !     fpar(7) == current rho (rhok = <r, r0>)
  !     fpar(8) == alpha
  !     fpar(9) == omega
  !
  !     Usage of the work space W
  !     w(:, 1) = r0, the initial residual vector
  !     w(:, 2) = r, current residual vector
  !     w(:, 3) = s
  !     w(:, 4) = t
  !     w(:, 5) = v
  !     w(:, 6) = p
  !     w(:, 7) = tmp, used in preconditioning, etc.
  !     w(:, 8) = delta x, the correction to the answer is accumulated
  !               here, so that the right-preconditioning may be applied
  !               at the end
  !-----------------------------------------------------------------------
  !     external routines used
  !
  real*8 ddot
  logical stopbis, brkdn
  external ddot, stopbis, brkdn
  !
  real*8 one
  parameter(one=1.0d0)
  !
  !     local variables
  !
  integer i
  real*8 alpha,beta,rho,omega
  logical lp, rp
  save lp, rp
  !
  !     where to go
  !
  if (ipar(1).gt.0) then
    !!goto (10, 20, 40, 50, 60, 70, 80, 90, 100, 110) ipar(10)
    SELECT CASE (ipar(10))
    CASE (1)
      GOTO 10
    CASE (2)
      GOTO 20
    CASE (3)
      GOTO 40
    CASE (4)
      GOTO 50
    CASE (5)
      GOTO 60
    CASE (6)
      GOTO 70
    CASE (7)
      GOTO 80
    CASE (8)
      GOTO 90
    CASE (9)
      GOTO 100
    CASE (10)
      GOTO 110
    END SELECT
  else if (ipar(1).lt.0) then
    goto 900
  endif
  !
  !     call the initialization routine
  !
  call bisinit(ipar,fpar,8*n,1,lp,rp,w)
  if (ipar(1).lt.0) return
  !
  !     perform a matvec to compute the initial residual
  !
  ipar(1) = 1
  ipar(8) = 1
  ipar(9) = 1 + n
  do i = 1, n
    w(i,1) = sol(i)
  enddo
  ipar(10) = 1
  return
10 ipar(7) = ipar(7) + 1
  ipar(13) = ipar(13) + 1
  do i = 1, n
    w(i,1) = rhs(i) - w(i,2)
  enddo
  fpar(11) = fpar(11) + n
  if (lp) then
    ipar(1) = 3
    ipar(10) = 2
    return
  endif
  !
20 if (lp) then
    do i = 1, n
      w(i,1) = w(i,2)
      w(i,6) = w(i,2)
    enddo
  else
    do i = 1, n
      w(i,2) = w(i,1)
      w(i,6) = w(i,1)
    enddo
  endif
  !
  fpar(7) = ddot(n,w,w)
  fpar(11) = fpar(11) + 2 * n
  fpar(5) = sqrt(fpar(7))
  fpar(3) = fpar(5)
  if (abs(ipar(3)).eq.2) then
    fpar(4) = fpar(1) * sqrt(ddot(n,rhs,rhs)) + fpar(2)
    fpar(11) = fpar(11) + 2 * n
  else if (ipar(3).ne.999) then
    fpar(4) = fpar(1) * fpar(3) + fpar(2)
  endif
  if (ipar(3).ge.0) fpar(6) = fpar(5)
  if (ipar(3).ge.0 .and. fpar(5).le.fpar(4) .and. ipar(3).ne.999) then
    goto 900
  endif
  !
  !     beginning of the iterations
  !
  !     Step (1), v = A p
30 if (rp) then
    ipar(1) = 5
    ipar(8) = 5*n+1
    if (lp) then
      ipar(9) = 4*n + 1
    else
      ipar(9) = 6*n + 1
    endif
    ipar(10) = 3
    return
  endif
  !
40 ipar(1) = 1
  if (rp) then
    ipar(8) = ipar(9)
  else
    ipar(8) = 5*n+1
  endif
  if (lp) then
    ipar(9) = 6*n + 1
  else
    ipar(9) = 4*n + 1
  endif
  ipar(10) = 4
  return
50 if (lp) then
    ipar(1) = 3
    ipar(8) = ipar(9)
    ipar(9) = 4*n + 1
    ipar(10) = 5
    return
  endif
  !
60 ipar(7) = ipar(7) + 1
  !
  !     step (2)
  alpha = ddot(n,w(1,1),w(1,5))
  fpar(11) = fpar(11) + 2 * n
  if (brkdn(alpha, ipar)) goto 900
  alpha = fpar(7) / alpha
  fpar(8) = alpha
  !
  !     step (3)
  do i = 1, n
    w(i,3) = w(i,2) - alpha * w(i,5)
  enddo
  fpar(11) = fpar(11) + 2 * n
  !
  !     Step (4): the second matvec -- t = A s
  !
  if (rp) then
    ipar(1) = 5
    ipar(8) = n+n+1
    if (lp) then
      ipar(9) = ipar(8)+n
    else
      ipar(9) = 6*n + 1
    endif
    ipar(10) = 6
    return
  endif
  !
70 ipar(1) = 1
  if (rp) then
    ipar(8) = ipar(9)
  else
    ipar(8) = n+n+1
  endif
  if (lp) then
    ipar(9) = 6*n + 1
  else
    ipar(9) = 3*n + 1
  endif
  ipar(10) = 7
  return
80 if (lp) then
    ipar(1) = 3
    ipar(8) = ipar(9)
    ipar(9) = 3*n + 1
    ipar(10) = 8
    return
  endif
90 ipar(7) = ipar(7) + 1
  !
  !     step (5)
  omega = ddot(n,w(1,4),w(1,4))
  fpar(11) = fpar(11) + n + n
  if (brkdn(omega,ipar)) goto 900
  omega = ddot(n,w(1,4),w(1,3)) / omega
  fpar(11) = fpar(11) + n + n
  if (brkdn(omega,ipar)) goto 900
  fpar(9) = omega
  alpha = fpar(8)
  !
  !     step (6) and (7)
  do i = 1, n
    w(i,7) = alpha * w(i,6) + omega * w(i,3)
    w(i,8) = w(i,8) + w(i,7)
    w(i,2) = w(i,3) - omega * w(i,4)
  enddo
  fpar(11) = fpar(11) + 6 * n + 1
  !
  !     convergence test
  if (ipar(3).eq.999) then
    ipar(1) = 10
    ipar(8) = 7*n + 1
    ipar(9) = 6*n + 1
    ipar(10) = 9
    return
  endif
  if (stopbis(n,ipar,2,fpar,w(1,2),w(1,7),one))  goto 900
100 if (ipar(3).eq.999.and.ipar(11).eq.1) goto 900
  !
  !     step (8): computing new p and rho
  !
  rho = fpar(7)
  fpar(7) = ddot(n,w(1,2),w(1,1))
  omega = fpar(9)
  beta = fpar(7) * fpar(8) / (fpar(9) * rho)
  do i = 1, n
    w(i,6) = w(i,2) + beta * (w(i,6) - omega * w(i,5))
  enddo
  fpar(11) = fpar(11) + 6 * n + 3
  if (brkdn(fpar(7),ipar)) goto 900
  !
  !     end of an iteration
  !
  goto 30
  !
  !     some clean up job to do
  !
900 if (rp) then
    if (ipar(1).lt.0) ipar(12) = ipar(1)
    ipar(1) = 5
    ipar(8) = 7*n + 1
    ipar(9) = ipar(8) - n
    ipar(10) = 10
    return
  endif
 
110 if (rp) then
    call tidycg(n,ipar,fpar,sol,w(1,7))
  else
    call tidycg(n,ipar,fpar,sol,w(1,8))
  endif
  !
  return
  !-----end-of-bcgstab

References bisinit(), brkdn(), ddot(), stopbis(), and tidycg().

Referenced by w3profsmd::w3xypfsnimp().

bisinit()

subroutine bisinit	(	integer, dimension(16)	ipar,
		real*8, dimension(16)	fpar,
		integer	wksize,
		integer	dsc,
		logical	lp,
		logical	rp,
		real*8, dimension(*)	wk
	)

Definition at line 2593 of file w3profsmd.F90.

  implicit none
  integer i,ipar(16),wksize,dsc
  logical lp,rp
  real*8  fpar(16),wk(*)
  !-----------------------------------------------------------------------
  !     some common initializations for the iterative solvers
  !-----------------------------------------------------------------------
  real*8 zero, one
  parameter(zero=0.0d0, one=1.0d0)
  !
  !     ipar(1) = -2 inidcate that there are not enough space in the work
  !     array
  !
  if (ipar(4).lt.wksize) then
    ipar(1) = -2
    ipar(4) = wksize
    return
  endif
  !
  if (ipar(2).gt.2) then
    lp = .true.
    rp = .true.
  else if (ipar(2).eq.2) then
    lp = .false.
    rp = .true.
  else if (ipar(2).eq.1) then
    lp = .true.
    rp = .false.
  else
    lp = .false.
    rp = .false.
  endif
  if (ipar(3).eq.0) ipar(3) = dsc
  !     .. clear the ipar elements used
  ipar(7) = 0
  ipar(8) = 0
  ipar(9) = 0
  ipar(10) = 0
  ipar(11) = 0
  ipar(12) = 0
  ipar(13) = 0
  !
  !     fpar(1) must be between (0, 1), fpar(2) must be positive,
  !     fpar(1) and fpar(2) can NOT both be zero
  !     Normally return ipar(1) = -4 to indicate any of above error
  !
  if (fpar(1).lt.zero .or. fpar(1).ge.one .or. fpar(2).lt.zero .or. &
       &     (fpar(1).eq.zero .and. fpar(2).eq.zero)) then
    if (ipar(1).eq.0) then
      ipar(1) = -4
      return
    else
      fpar(1) = 1.0d-6
      fpar(2) = 1.0d-16
    endif
  endif
  !     .. clear the fpar elements
  do i = 3, 10
    fpar(i) = zero
  enddo
  if (fpar(11).lt.zero) fpar(11) = zero
  !     .. clear the used portion of the work array to zero
  do i = 1, wksize
    wk(i) = zero
  enddo
  !
  return
  !-----end-of-bisinit

Referenced by bcgstab().

brkdn()

logical function brkdn	(	real*8	alpha,
		integer, dimension(16)	ipar
	)

Definition at line 2557 of file w3profsmd.F90.

  implicit none
  integer ipar(16)
  real*8 alpha, beta, zero, one
  parameter(zero=0.0d0, one=1.0d0)
  !-----------------------------------------------------------------------
  !     test whether alpha is zero or an abnormal number, if yes,
  !     this routine will return .true.
  !
  !     If alpha == 0, ipar(1) = -3,
  !     if alpha is an abnormal number, ipar(1) = -9.
  !-----------------------------------------------------------------------
  brkdn = .false.
  if (alpha.gt.zero) then
    beta = one / alpha
    if (.not. beta.gt.zero) then
      brkdn = .true.
      ipar(1) = -9
    endif
  else if (alpha.lt.zero) then
    beta = one / alpha
    if (.not. beta.lt.zero) then
      brkdn = .true.
      ipar(1) = -9
    endif
  else if (alpha.eq.zero) then
    brkdn = .true.
    ipar(1) = -3
  else
    brkdn = .true.
    ipar(1) = -9
  endif
  return

Referenced by bcgstab().

daxpy()

subroutine daxpy	(	integer	n,
		double precision	da,
		double precision, dimension(1)	dx,
		integer	incx,
		double precision, dimension(1)	dy,
		integer	incy
	)

Definition at line 4641 of file w3profsmd.F90.

  !
  !     constant times a vector plus a vector.
  !     uses unrolled loops for increments equal to one.
  !     jack dongarra, linpack, 3/11/78.
  !
  double precision dx(1),dy(1),da
  integer i,incx,incy,ix,iy,m,mp1,n
  !
  if(n.le.0)return
  if (abs(da) .lt. tiny(1.d0)) 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  i = 1,n
    dy(iy) = dy(iy) + da*dx(ix)
    ix = ix + incx
    iy = iy + incy
  end do
  return
  !
  !        code for both increments equal to 1
  !
  !
  !        clean-up loop
  !
20 m = mod(n,4)
  if( m .eq. 0 ) go to 40
  do i = 1,m
    dy(i) = dy(i) + da*dx(i)
  end do
  if( n .lt. 4 ) return
40 mp1 = m + 1
  do i = mp1,n,4
    dy(i) = dy(i) + da*dx(i)
    dy(i + 1) = dy(i + 1) + da*dx(i + 1)
    dy(i + 2) = dy(i + 2) + da*dx(i + 2)
    dy(i + 3) = dy(i + 3) + da*dx(i + 3)
  end do
  return

Referenced by pgmres().

ddot()

double precision function ddot	(	integer	n,
		double precision, dimension(*)	dx,
		double precision, dimension(*)	dy
	)

Definition at line 4613 of file w3profsmd.F90.

  !
  !     forms the dot product of two vectors.
  !     uses unrolled loops for increments equal to one.
  !     jack dongarra, linpack, 3/11/78.
  !
  double precision dx(*),dy(*),dtemp
  integer i,m,mp1,n
  !
  ddot = 0.0d0
  dtemp = 0.0d0
  if(n.le.0)return
 
20 m = mod(n,5)
  if( m .eq. 0 ) go to 40
  do i = 1,m
    dtemp = dtemp + dx(i)*dy(i)
  end do
  if( n .lt. 5 ) go to 60
40 mp1 = m + 1
  do i = mp1,n,5
    dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + &
         &   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
  end do
60 ddot = dtemp
  return

Referenced by bcgstab(), mgsro(), pgmres(), and stopbis().

dlassq()

subroutine dlassq	(	integer	N,
		double precision, dimension( * )	X,
		double precision	SCALE,
		double precision	SUMSQ
	)

Definition at line 4568 of file w3profsmd.F90.

  !
  ! -- LAPACK auxiliary routine (version 3.1) --
  ! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  ! November 2006
  INTEGER N
  DOUBLE PRECISION SCALE, SUMSQ
  DOUBLE PRECISION X( * )
  !
  ! DLASSQ returns the values scl and smsq such that
  !
  ! ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
  !
  ! where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is
  ! assumed to be non-negative and scl returns the value
  !
  ! scl = max( scale, abs( x( i ) ) ).
  !
  ! SCALE (input/output) DOUBLE PRECISION
  ! On entry, the value scale in the equation above.
  ! On exit, SCALE is overwritten with scl , the scaling factor
  ! for the sum of squares.
  DOUBLE PRECISION ZERO
  parameter( zero = 0.0d+0 )
  INTEGER IX
  DOUBLE PRECISION ABSXI
  INTRINSIC abs
  !
  IF( n.GT.0 ) THEN
    DO ix = 1, 1 + ( n-1 )
      IF( x( ix ).NE.zero ) THEN
        absxi = abs( x( ix ) )
        IF( scale.LT.absxi ) THEN
          sumsq = 1 + sumsq*( scale / absxi )**2
          scale = absxi
        ELSE
          sumsq = sumsq + ( absxi / scale )**2
        END IF
      END IF
    END DO
  END IF
  RETURN

dnrm2()

double precision function dnrm2	(	integer	N,
		double precision, dimension(*)	X
	)

Definition at line 4499 of file w3profsmd.F90.

  !     .. Scalar Arguments ..
  INTEGER N
  !     ..
  !     .. Array Arguments ..
  DOUBLE PRECISION X(*)
  !     ..
  !
  !  Purpose
  !  =======
  !
  !  DNRM2 returns the euclidean norm of a vector via the function
  !  name, so that
  !
  !     DNRM2 := sqrt( x'*x )
  !
  !  Further Details
  !  ===============
  !
  !  -- This version written on 25-October-1982.
  !     Modified on 14-October-1993 to inline the call to DLASSQ.
  !     Sven Hammarling, Nag Ltd.
  !
  !  =====================================================================
  !
  !     .. Parameters ..
  DOUBLE PRECISION ONE,ZERO
  parameter(one=1.0d+0,zero=0.0d+0)
  !     ..
  !     .. Local Scalars ..
  DOUBLE PRECISION ABSXI,NORM,SCALE,SSQ
  INTEGER IX
  !     ..
  !     .. Intrinsic Functions ..
  INTRINSIC abs,sqrt
  !     ..
  IF (n.LT.1 ) THEN
    norm = zero
  ELSE IF (n.EQ.1) THEN
    norm = abs(x(1))
  ELSE
    scale = zero
    ssq = one
    !        The following loop is equivalent to this call to the LAPACK
    !        auxiliary routine:
    !        CALL DLASSQ( N, X, SCALE, SSQ )
    !
    DO ix = 1,1 + (n-1)
      IF (x(ix).NE.zero) THEN
        absxi = abs(x(ix))
        IF (scale.LT.absxi) THEN
          ssq = one + ssq* (scale/absxi)**2
          scale = absxi
        ELSE
          ssq = ssq + (absxi/scale)**2
        END IF
      END IF
    end do
    norm = scale*sqrt(ssq)
  END IF
  !
  dnrm2 = norm
  RETURN
  !
  !     End of DNRM2.
  !

Referenced by pgmres().

givens()

subroutine givens	(	real*8	x,
		real*8	y,
		real*8	c,
		real*8	s
	)

Definition at line 2423 of file w3profsmd.F90.

  implicit none
 
  real*8   :: x,y,c,s
  !-----------------------------------------------------------------------
  !     Given x and y, this subroutine generates a Givens' rotation c, s.
  !     And apply the rotation on (x,y) ==> (sqrt(x**2 + y**2), 0).
  !     (See P 202 of "matrix computation" by Golub and van Loan.)
  !-----------------------------------------------------------------------
  real*8   :: t,one,zero
  parameter(zero=0.0d0,one=1.0d0)
  !
  if (x.eq.zero .and. y.eq.zero) then
    c = one
    s = zero
  else if (abs(y).gt.abs(x)) then
    t = x / y
    x = sqrt(one+t*t)
    s = sign(one / x, y)
    c = t*s
  else if (abs(y).le.abs(x)) then
    t = y / x
    y = sqrt(one+t*t)
    c = sign(one / y, x)
    s = t*c
  else
    !
    !     X or Y must be an invalid floating-point number, set both to zero
    !
    x = zero
    y = zero
    c = one
    s = zero
  endif
  x = abs(x*y)
  !
  !     end of givens
  !
  return

ilu0()

subroutine ilu0	(	integer	n,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(*)	ia,
		real*8, dimension(*)	alu,
		integer, dimension(*)	jlu,
		integer, dimension(*)	ju,
		integer, dimension(n)	iw,
		integer	ierr
	)

Definition at line 4177 of file w3profsmd.F90.

 
  !implicit real*8 (a-h,o-z)
  real*8 a(*), alu(*), tl
  integer n, ju0, ii, jj, i, j, jcol, js, jf, jm, jrow, jw, ierr
  integer ja(*), ia(*), ju(*), jlu(*), iw(n)
  !
  !-----------------------------------------------------------------------
  ju0 = n+2
  jlu(1) = ju0 !!!
  iw = 0
  do ii = 1, n
    js = ju0
    do j=ia(ii),ia(ii+1)-1
      jcol = ja(j)
      if (jcol .eq. ii) then
        alu(ii) = a(j)
        iw(jcol) = ii
        ju(ii)  = ju0 !!!
      else
        alu(ju0) = a(j)
        jlu(ju0) = ja(j)
        iw(jcol) = ju0
        ju0 = ju0+1
      endif
    end do
    jlu(ii+1) = ju0 !!!
    jf = ju0-1
    jm = ju(ii)-1
    !     exit if diagonal element is reached.
    do j=js, jm
      jrow = jlu(j)
      tl = alu(j)*alu(jrow)
      alu(j) = tl
      !     perform  linear combination
      do jj = ju(jrow), jlu(jrow+1)-1
        jw = iw(jlu(jj))
        if (jw .ne. 0) then
          alu(jw) = alu(jw) - tl*alu(jj)
          !                   write(*,*) ii, jw, jj
        end if
      end do
    end do
    !     invert  and store diagonal element.
    if (abs(alu(ii)) .lt. tiny(1.)) goto 600
    alu(ii) = 1.0d0/alu(ii)
    !     reset pointer iw to zero
    iw(ii) = 0
    do i = js, jf
      iw(jlu(i)) = 0
    end do
  end do
 
  ierr = 0
  return
 
  !     zero pivot :
600 ierr = ii
  return

Referenced by w3profsmd::w3xypfsnimp().

ilut()

subroutine ilut	(	integer	n,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(n+1)	ia,
		integer	lfil,
		real*8	droptol,
		real*8, dimension(*)	alu,
		integer, dimension(*)	jlu,
		integer, dimension(n)	ju,
		integer	iwk,
		real*8, dimension(n+1)	w,
		integer, dimension(2*n)	jw,
		integer	ierr
	)

Definition at line 3827 of file w3profsmd.F90.

  !-----------------------------------------------------------------------
  implicit none
  integer n
  real*8 a(*),alu(*),w(n+1),droptol
  integer ja(*),ia(n+1),jlu(*),ju(n),jw(2*n),lfil,iwk,ierr
  !----------------------------------------------------------------------*
  !                      *** ILUT preconditioner ***                     *
  !      incomplete LU factorization with dual truncation mechanism      *
  !----------------------------------------------------------------------*
  !     Author: Yousef Saad *May, 5, 1990, Latest revision, August 1996  *
  !----------------------------------------------------------------------*
  ! PARAMETERS
  !-----------
  !
  ! on entry:
  !==========
  ! n       = integer. The row dimension of the matrix A. The matrix
  !
  ! a,ja,ia = matrix stored in Compressed Sparse Row format.
  !
  ! lfil    = integer. The fill-in parameter. Each row of L and each row
  !           of U will have a maximum of lfil elements (excluding the
  !           diagonal element). lfil must be .ge. 0.
  !           ** WARNING: THE MEANING OF LFIL HAS CHANGED WITH RESPECT TO
  !           EARLIER VERSIONS.
  !
  ! droptol = real*8. Sets the threshold for dropping small terms in the
  !           factorization. See below for details on dropping strategy.
  !
  !
  ! iwk     = integer. The lengths of arrays alu and jlu. If the arrays
  !           are not big enough to store the ILU factorizations, ilut
  !           will stop with an error message.
  !
  ! On return:
  !===========
  !
  ! alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
  !           the L and U factors together. The diagonal (stored in
  !           alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
  !           contains the i-th row of L (excluding the diagonal entry=1)
  !           followed by the i-th row of U.
  !
  ! ju      = integer array of length n containing the pointers to
  !           the beginning of each row of U in the matrix alu,jlu.
  !
  ! ierr    = integer. Error message with the following meaning.
  !           ierr  = 0    --> successful return.
  !           ierr .gt. 0  --> zero pivot encountered at step number ierr.
  !           ierr  = -1   --> Error. input matrix may be wrong.
  !                            (The elimination process has generated a
  !                            row in L or U whose length is .gt.  n.)
  !           ierr  = -2   --> The matrix L overflows the array al.
  !           ierr  = -3   --> The matrix U overflows the array alu.
  !           ierr  = -4   --> Illegal value for lfil.
  !           ierr  = -5   --> zero row encountered.
  !
  ! work arrays:
  !=============
  ! jw      = integer work array of length 2*n.
  ! w       = real work array of length n+1.
  !
  !----------------------------------------------------------------------
  ! w, ju (1:n) store the working array [1:ii-1 = L-part, ii:n = u]
  ! jw(n+1:2n)  stores nonzero indicators
  !
  ! Notes:
  ! ------
  ! The diagonal elements of the input matrix must be  nonzero (at least
  ! 'structurally').
  !
  !----------------------------------------------------------------------*
  !---- Dual drop strategy works as follows.                             *
  !                                                                      *
  !     1) Theresholding in L and U as set by droptol. Any element whose *
  !        magnitude is less than some tolerance (relative to the abs    *
  !        value of diagonal element in u) is dropped.                   *
  !                                                                      *
  !     2) Keeping only the largest lfil elements in the i-th row of L   *
  !        and the largest lfil elements in the i-th row of U (excluding *
  !        diagonal elements).                                           *
  !                                                                      *
  ! Flexibility: one  can use  droptol=0  to get  a strategy  based on   *
  ! keeping  the largest  elements in  each row  of L  and U.   Taking   *
  ! droptol .ne.  0 but lfil=n will give  the usual threshold strategy   *
  ! (however, fill-in is then mpredictible).                             *
  !----------------------------------------------------------------------*
  !     locals
  integer ju0,k,j1,j2,j,ii,i,lenl,lenu,jj,jrow,jpos,lenn
  real*8 tnorm, t, abs, s, fact
  if (lfil .lt. 0) goto 998
  !-----------------------------------------------------------------------
  !     initialize ju0 (points to next element to be added to alu,jlu)
  !     and pointer array.
  !-----------------------------------------------------------------------
  ju0 = n+2
  jlu(1) = ju0
  !
  !     initialize nonzero indicator array.
  !
  do  j=1,n
    jw(n+j)  = 0
  end do
  !-----------------------------------------------------------------------
  !     beginning of main loop.
  !-----------------------------------------------------------------------
  do ii = 1, n
 
    j1 = ia(ii)
    j2 = ia(ii+1) - 1
 
    tnorm = 0.0d0
    do k=j1,j2
      tnorm = tnorm+abs(a(k))
    end do
    if (abs(tnorm) .lt. tiny(1.)) goto 999
 
    tnorm = tnorm/real(j2-j1+1)
    !
    !     unpack L-part and U-part of row of A in arrays w
    !
    lenu = 1
    lenl = 0
    jw(ii) = ii
    w(ii) = 0.0
    jw(n+ii) = ii
    !
    do j = j1, j2
      k = ja(j)
      t = a(j)
      if (k .lt. ii) then
        lenl = lenl+1
        jw(lenl) = k
        w(lenl) = t
        jw(n+k) = lenl
      else if (k .eq. ii) then
        w(ii) = t
      else
        lenu = lenu+1
        jpos = ii+lenu-1
        jw(jpos) = k
        w(jpos) = t
        jw(n+k) = jpos
      endif
    end do
    jj = 0
    lenn = 0
    !
    !     eliminate previous rows
    !
150 jj = jj+1
    if (jj .gt. lenl) goto 160
    !-----------------------------------------------------------------------
    !     in order to do the elimination in the correct order we must select
    !     the smallest column index among jw(k), k=jj+1, ..., lenl.
    !-----------------------------------------------------------------------
    jrow = jw(jj)
    k = jj
    !
    !     determine smallest column index
    !
    do j=jj+1,lenl
      if (jw(j) .lt. jrow) then
        jrow = jw(j)
        k = j
      endif
    end do
    !
    if (k .ne. jj) then
      !     exchange in jw
      j = jw(jj)
      jw(jj) = jw(k)
      jw(k) = j
      !     exchange in jr
      jw(n+jrow) = jj
      jw(n+j) = k
      !     exchange in w
      s = w(jj)
      w(jj) = w(k)
      w(k) = s
    endif
    !
    !     zero out element in row by setting jw(n+jrow) to zero.
    !
    jw(n+jrow) = 0
    !
    !     get the multiplier for row to be eliminated (jrow).
    !
    fact = w(jj)*alu(jrow)
    if (abs(fact) .le. droptol) goto 150
    !
    !     combine current row and row jrow
    !
    do  k = ju(jrow), jlu(jrow+1)-1
      s = fact*alu(k)
      j = jlu(k)
      jpos = jw(n+j)
      if (j .ge. ii) then
        !
        !     dealing with upper part.
        !
        if (jpos .eq. 0) then
          !
          !     this is a fill-in element
          !
          lenu = lenu+1
          if (lenu .gt. n) goto 995
          i = ii+lenu-1
          jw(i) = j
          jw(n+j) = i
          w(i) = - s
        else
          !
          !     this is not a fill-in element
          !
          w(jpos) = w(jpos) - s
 
        endif
      else
        !
        !     dealing  with lower part.
        !
        if (jpos .eq. 0) then
          !
          !     this is a fill-in element
          !
          lenl = lenl+1
          if (lenl .gt. n) goto 995
          jw(lenl) = j
          jw(n+j) = lenl
          w(lenl) = - s
        else
          !
          !     this is not a fill-in element
          !
          w(jpos) = w(jpos) - s
        endif
      endif
    end do
    !
    !     store this pivot element -- (from left to right -- no danger of
    !     overlap with the working elements in L (pivots).
    !
    lenn = lenn+1
    w(lenn) = fact
    jw(lenn)  = jrow
    goto 150
160 continue
    !
    !     reset double-pointer to zero (U-part)
    !
    do k=1, lenu
      jw(n+jw(ii+k-1)) = 0
    end do
    !
    !     update L-matrix
    !
    lenl = lenn
    lenn = min0(lenl,lfil)
    !
    !     sort by quick-split
    !
    call qsplit (w,jw,lenl,lenn)
    !
    !     store L-part
    !
    do k=1, lenn
      if (ju0 .gt. iwk) goto 996
      alu(ju0) =  w(k)
      jlu(ju0) =  jw(k)
      ju0 = ju0+1
    end do
    !
    !     save pointer to beginning of row ii of U
    !
    ju(ii) = ju0
    !
    !     update U-matrix -- first apply dropping strategy
    !
    lenn = 0
    do k=1, lenu-1
      if (abs(w(ii+k)) .gt. droptol*tnorm) then
        lenn = lenn+1
        w(ii+lenn) = w(ii+k)
        jw(ii+lenn) = jw(ii+k)
      endif
    enddo
    lenu = lenn+1
    lenn = min0(lenu,lfil)
    !
    call qsplit (w(ii+1), jw(ii+1), lenu-1,lenn)
    !
    !     copy
    !
    t = abs(w(ii))
    if (lenn + ju0 .gt. iwk) goto 997
    do k=ii+1,ii+lenn-1
      jlu(ju0) = jw(k)
      alu(ju0) = w(k)
      t = t + abs(w(k) )
      ju0 = ju0+1
    end do
    !
    !     store inverse of diagonal element of u
    !
    !2do check if it works ... after correction ...
    if (abs(w(ii)) .lt.  tiny(1.d0)) w(ii) = (0.0001d0 + droptol)*tnorm
    !
    alu(ii) = 1.0d0/ w(ii)
    !
    !     update pointer to beginning of next row of U.
    !
    jlu(ii+1) = ju0
    !-----------------------------------------------------------------------
    !     end main loop
    !-----------------------------------------------------------------------
  end do
  ierr = 0
  return
  !
  !     incomprehensible error. Matrix must be wrong.
  !
995 ierr = -1
  return
  !
  !     insufficient storage in L.
  !
996 ierr = -2
  return
  !
  !     insufficient storage in U.
  !
997 ierr = -3
  return
  !
  !     illegal lfil entered.
  !
998 ierr = -4
  return
  !
  !     zero row encountered
  !
999 ierr = -5
  return
  !----------------end-of-ilut--------------------------------------------
  !-----------------------------------------------------------------------

References qsplit().

implu()

subroutine implu	(	integer	np,
		real*8	umm,
		real*8	beta,
		real*8, dimension(*)	ypiv,
		real*8, dimension(*)	u,
		logical, dimension(*)	permut,
		logical	full
	)

Definition at line 2342 of file w3profsmd.F90.

  real*8 umm,beta,ypiv(*),u(*),x, xpiv
  logical full, perm, permut(*)
  integer np,k,npm1
  !-----------------------------------------------------------------------
  !     performs implicitly one step of the lu factorization of a
  !     banded hessenberg matrix.
  !-----------------------------------------------------------------------
  if (np .le. 1) goto 12
  npm1 = np - 1
  !
  !     -- perform  previous step of the factorization-
  !
  do k=1,npm1
    if (.not. permut(k)) goto 5
    x=u(k)
    u(k) = u(k+1)
    u(k+1) = x
5   u(k+1) = u(k+1) - ypiv(k)*u(k)
  end do
  !-----------------------------------------------------------------------
  !     now determine pivotal information to be used in the next call
  !-----------------------------------------------------------------------
12 umm = u(np)
  perm = (beta .gt. abs(umm))
  if (.not. perm) goto 4
  xpiv = umm / beta
  u(np) = beta
  goto 8
4 xpiv = beta/umm
8 permut(np) = perm
  ypiv(np) = xpiv
  if (.not. full) return
  !     shift everything up if full...
  do k=1,npm1
    ypiv(k) = ypiv(k+1)
    permut(k) = permut(k+1)
  end do
  return
  !-----end-of-implu

ldsol()

subroutine ldsol	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	al,
		integer, dimension(*)	jal
	)

Definition at line 3258 of file w3profsmd.F90.

  integer n, jal(*)
  real*8 x(n), y(n), al(*)
  !-----------------------------------------------------------------------
  !     Solves L x = y    L = triangular. MSR format
  !-----------------------------------------------------------------------
  ! solves a (non-unit) lower triangular system by standard (sequential)
  ! forward elimination - matrix stored in MSR format
  ! with diagonal elements already inverted (otherwise do inversion,
  ! al(1:n) = 1.0/al(1:n),  before calling ldsol).
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real array containg the right hand side.
  !
  ! al,
  ! jal,   = Lower triangular matrix stored in Modified Sparse Row
  !          format.
  !
  ! On return:
  !-----------
  !        x = The solution of  L x = y .
  !--------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8 t
  !-----------------------------------------------------------------------
  x(1) = y(1)*al(1)
  do  k = 2, n
    t = y(k)
    do  j = jal(k), jal(k+1)-1
      t = t - al(j)*x(jal(j))
    end do
    x(k) = al(k)*t
  end do
  return
  !----------end-of-ldsol-------------------------------------------------
  !-----------------------------------------------------------------------

ldsolc()

subroutine ldsolc	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	al,
		integer, dimension(*)	jal
	)

Definition at line 3345 of file w3profsmd.F90.

  integer n, jal(*)
  real*8 x(n), y(n), al(*)
  !-----------------------------------------------------------------------
  !    Solves     L x = y ;    L = nonunit Low. Triang. MSC format
  !-----------------------------------------------------------------------
  ! solves a (non-unit) lower triangular system by standard (sequential)
  ! forward elimination - matrix stored in Modified Sparse Column format
  ! with diagonal elements already inverted (otherwise do inversion,
  ! al(1:n) = 1.0/al(1:n),  before calling ldsol).
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real array containg the right hand side.
  !
  ! al,
  ! jal,
  ! ial,    = Lower triangular matrix stored in Modified Sparse Column
  !           format.
  !
  ! On return:
  !-----------
  !         x = The solution of  L x = y .
  !--------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8 t
  !-----------------------------------------------------------------------
  do k=1,n
    x(k) = y(k)
  end do
  do k = 1, n
    x(k) = x(k)*al(k)
    t = x(k)
    do j = jal(k), jal(k+1)-1
      x(jal(j)) = x(jal(j)) - t*al(j)
    end do
  end do
  !
  return
  !----------end-of-lsolc------------------------------------------------
  !-----------------------------------------------------------------------

ldsoll()

subroutine ldsoll	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	al,
		integer, dimension(*)	jal,
		integer	nlev,
		integer, dimension(n)	lev,
		integer, dimension(nlev+1)	ilev
	)

Definition at line 3392 of file w3profsmd.F90.

  integer n, nlev, jal(*), ilev(nlev+1), lev(n)
  real*8 x(n), y(n), al(*)
  !-----------------------------------------------------------------------
  !    Solves L x = y    L = triangular. Uses LEVEL SCHEDULING/MSR format
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real array containg the right hand side.
  !
  ! al,
  ! jal,   = Lower triangular matrix stored in Modified Sparse Row
  !          format.
  ! nlev   = number of levels in matrix
  ! lev    = integer array of length n, containing the permutation
  !          that defines the levels in the level scheduling ordering.
  ! ilev   = pointer to beginning of levels in lev.
  !          the numbers lev(i) to lev(i+1)-1 contain the row numbers
  !          that belong to level number i, in the level shcheduling
  !          ordering.
  !
  ! On return:
  !-----------
  !        x = The solution of  L x = y .
  !--------------------------------------------------------------------
  integer ii, jrow, i, k
  real*8 t
  !
  !     outer loop goes through the levels. (SEQUENTIAL loop)
  !
  do ii=1, nlev
    !
    !     next loop executes within the same level. PARALLEL loop
    !
    do i=ilev(ii), ilev(ii+1)-1
      jrow = lev(i)
      !
      ! compute inner product of row jrow with x
      !
      t = y(jrow)
      do k=jal(jrow), jal(jrow+1)-1
        t = t - al(k)*x(jal(k))
      end do
      x(jrow) = t*al(jrow)
    end do
  end do
  return
  !-----------------------------------------------------------------------

lsol()

subroutine lsol	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	al,
		integer, dimension(*)	jal,
		integer, dimension(n+1)	ial
	)

Definition at line 3215 of file w3profsmd.F90.

  integer n, jal(*),ial(n+1)
  real*8  x(n), y(n), al(*)
  !-----------------------------------------------------------------------
  !   solves    L x = y ; L = lower unit triang. /  CSR format
  !-----------------------------------------------------------------------
  ! solves a unit lower triangular system by standard (sequential )
  ! forward elimination - matrix stored in CSR format.
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real array containg the right side.
  !
  ! al,
  ! jal,
  ! ial,    = Lower triangular matrix stored in compressed sparse row
  !          format.
  !
  ! On return:
  !-----------
  !          x  = The solution of  L x  = y.
  !--------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8  t
  !-----------------------------------------------------------------------
  x(1) = y(1)
  do  k = 2, n
    t = y(k)
    do j = ial(k), ial(k+1)-1
      t = t-al(j)*x(jal(j))
    end do
    x(k) = t
  end do
  !
  return
  !----------end-of-lsol--------------------------------------------------
  !-----------------------------------------------------------------------

lsolc()

subroutine lsolc	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	al,
		integer, dimension(*)	jal,
		integer, dimension(*)	ial
	)

Definition at line 3301 of file w3profsmd.F90.

  integer n, jal(*),ial(*)
  real*8  x(n), y(n), al(*)
  !-----------------------------------------------------------------------
  !       SOLVES     L x = y ;    where L = unit lower trang. CSC format
  !-----------------------------------------------------------------------
  ! solves a unit lower triangular system by standard (sequential )
  ! forward elimination - matrix stored in CSC format.
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real*8 array containg the right side.
  !
  ! al,
  ! jal,
  ! ial,    = Lower triangular matrix stored in compressed sparse column
  !          format.
  !
  ! On return:
  !-----------
  !         x  = The solution of  L x  = y.
  !-----------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8 t
  !-----------------------------------------------------------------------
  do k=1,n
    x(k) = y(k)
  end do
  do k = 1, n-1
    t = x(k)
    do j = ial(k), ial(k+1)-1
      x(jal(j)) = x(jal(j)) - t*al(j)
    end do
  end do
  !
  return
  !----------end-of-lsolc-------------------------------------------------
  !-----------------------------------------------------------------------

lusol()

subroutine lusol	(	integer	n,
		real*8, dimension(n)	y,
		real*8, dimension(n)	x,
		real*8, dimension(*)	alu,
		integer, dimension(*)	jlu,
		integer, dimension(*)	ju
	)

Definition at line 3621 of file w3profsmd.F90.

  implicit none
 
  integer    :: n, jlu(*), ju(*)
  real*8     :: x(n), y(n), alu(*)
 
  !-----------------------------------------------------------------------
  integer    :: i,k
  !
  ! forward solve
  !
  do i = 1, n
    x(i) = y(i)
    do k=jlu(i),ju(i)-1
      x(i) = x(i) - alu(k)* x(jlu(k))
    end do
  end do
  do i = n, 1, -1
    do k=ju(i),jlu(i+1)-1
      x(i) = x(i) - alu(k)*x(jlu(k))
    end do
    x(i) = alu(i)*x(i)
  end do
  !
  return
  !----------------end of lusol ------------------------------------------

Referenced by pgmres(), and runrc().

lutsol()

subroutine lutsol	(	integer	n,
		real*8, dimension(n)	y,
		real*8, dimension(n)	x,
		real*8, dimension(*)	alu,
		integer, dimension(*)	jlu,
		integer, dimension(*)	ju
	)

Definition at line 3649 of file w3profsmd.F90.

  implicit none
 
  integer     :: n, jlu(*), ju(*)
  real*8      :: x(n), y(n), alu(*)
 
  !-----------------------------------------------------------------------
  ! local variables
  !
  integer      :: i,k
  !
  do i = 1, n
    x(i) = y(i)
  end do
  !
  ! forward solve (with U^T)
  !
  do i = 1, n
    x(i) = x(i) * alu(i)
    do k=ju(i),jlu(i+1)-1
      x(jlu(k)) = x(jlu(k)) - alu(k)* x(i)
    end do
  end do
  !
  !     backward solve (with L^T)
  !
  do i = n, 1, -1
    do k=jlu(i),ju(i)-1
      x(jlu(k)) = x(jlu(k)) - alu(k)*x(i)
    end do
  end do
  !
  return
  !----------------end of lutsol -----------------------------------------
  !-----------------------------------------------------------------------

Referenced by runrc().

mgsro()

subroutine mgsro	(	logical	full,
		integer	lda,
		integer	n,
		integer	m,
		integer	ind,
		real*8	ops,
		real*8, dimension(lda,m)	vec,
		real*8, dimension(m)	hh,
		integer	ierr
	)

Definition at line 2664 of file w3profsmd.F90.

  implicit none
  logical full
  integer lda,m,n,ind,ierr
  real*8  ops,hh(m),vec(lda,m)
  !-----------------------------------------------------------------------
  !     MGSRO  -- Modified Gram-Schmidt procedure with Selective Re-
  !               Orthogonalization
  !     The ind'th vector of VEC is orthogonalized against the rest of
  !     the vectors.
  !
  !     The test for performing re-orthogonalization is performed for
  !     each indivadual vectors. If the cosine between the two vectors
  !     is greater than 0.99 (REORTH = 0.99**2), re-orthogonalization is
  !     performed. The norm of the 'new' vector is kept in variable NRM0,
  !     and updated after operating with each vector.
  !
  !     full   -- .ture. if it is necessary to orthogonalize the ind'th
  !               against all the vectors vec(:,1:ind-1), vec(:,ind+2:m)
  !               .false. only orthogonalize againt vec(:,1:ind-1)
  !     lda    -- the leading dimension of VEC
  !     n      -- length of the vector in VEC
  !     m      -- number of vectors can be stored in VEC
  !     ind    -- index to the vector to be changed
  !     ops    -- operation counts
  !     vec    -- vector of LDA X M storing the vectors
  !     hh     -- coefficient of the orthogonalization
  !     ierr   -- error code
  !               0 : successful return
  !               -1: zero input vector
  !               -2: input vector contains abnormal numbers
  !               -3: input vector is a linear combination of others
  !
  !     External routines used: real*8 ddot
  !-----------------------------------------------------------------------
  integer i,k
  real*8  nrm0, nrm1, fct, thr, ddot, zero, one, reorth
  parameter(zero=0.0d0, one=1.0d0, reorth=0.98d0)
  external ddot
  !
  !     compute the norm of the input vector
  !
  nrm0 = ddot(n,vec(1,ind),vec(1,ind))
  ops = ops + n + n
  thr = nrm0 * reorth
  if (nrm0.le.zero) then
    ierr = - 1
    return
  else if (nrm0.gt.zero .and. one/nrm0.gt.zero) then
    ierr = 0
  else
    ierr = -2
    return
  endif
  !
  !     Modified Gram-Schmidt loop
  !
  if (full) then
    do i = ind+1, m
      fct = ddot(n,vec(1,ind),vec(1,i))
      hh(i) = fct
      do k = 1, n
        vec(k,ind) = vec(k,ind) - fct * vec(k,i)
      end do
      ops = ops + 4 * n + 2
      if (fct*fct.gt.thr) then
        fct = ddot(n,vec(1,ind),vec(1,i))
        hh(i) = hh(i) + fct
        do k = 1, n
          vec(k,ind) = vec(k,ind) - fct * vec(k,i)
        end do
        ops = ops + 4*n + 1
      endif
      nrm0 = nrm0 - hh(i) * hh(i)
      if (nrm0.lt.zero) nrm0 = zero
      thr = nrm0 * reorth
    end do
  endif
  !
  do i = 1, ind-1
    fct = ddot(n,vec(1,ind),vec(1,i))
    hh(i) = fct
    do k = 1, n
      vec(k,ind) = vec(k,ind) - fct * vec(k,i)
    end do
    ops = ops + 4 * n + 2
    if (fct*fct.gt.thr) then
      fct = ddot(n,vec(1,ind),vec(1,i))
      hh(i) = hh(i) + fct
      do k = 1, n
        vec(k,ind) = vec(k,ind) - fct * vec(k,i)
      end do
      ops = ops + 4*n + 1
    endif
    nrm0 = nrm0 - hh(i) * hh(i)
    if (nrm0.lt.zero) nrm0 = zero
    thr = nrm0 * reorth
  end do
  !
  !     test the resulting vector
  !
  nrm1 = sqrt(ddot(n,vec(1,ind),vec(1,ind)))
  ops = ops + n + n
  hh(ind) = nrm1    ! statement label 75
  if (nrm1.le.zero) then
    ierr = -3
    return
  endif
  !
  !     scale the resulting vector
  !
  fct = one / nrm1
  do k = 1, n
    vec(k,ind) = vec(k,ind) * fct
  end do
  ops = ops + n + 1
  !
  !     normal return
  !
  ierr = 0
  return
  !     end subroutine mgsro

References ddot().

pgmres()

subroutine pgmres	(	integer	n,
		integer	im,
		real*8, dimension(*)	rhs,
		real*8, dimension(*)	sol,
		real*8	eps,
		integer	maxits,
		real*8, dimension(nnz)	aspar,
		integer	nnz,
		integer, dimension(n+1)	ia,
		integer, dimension(nnz)	ja,
		real*8, dimension(nnz+1)	alu,
		integer, dimension(nnz+1)	jlu,
		integer, dimension(n)	ju,
		real*8, dimension(n,im+1)	vv,
		integer	ierr
	)

Definition at line 4240 of file w3profsmd.F90.

  !-----------------------------------------------------------------------
  !       use datapool, only : nnz, ia, ja, alu, jlu, ju, vv, aspar!, rhs, sol
  implicit none
 
  integer :: n, im, maxits, ierr, nnz
 
  integer :: ja(nnz), ia(n+1)
  integer :: jlu(nnz+1), ju(n)
  real*8  :: vv(n,im+1), alu(nnz+1)
  real*8  :: aspar(nnz)
 
  real*8  :: rhs(*), sol(*)
 
  real*8  :: eps
  real*8  :: eps1, epsmac, gam, t, ddot, dnrm2, ro, tl
 
  integer :: i,i1,j,jj,k,k1,iii,ii,ju0
  integer :: its,jrow,jcol,jf,jm,js,jw
 
  real*8  :: hh(im+1,im), c(im), s(im), rs(im+1)
  real*8  :: iw(n)
 
  logical :: lblas = .false.      ! use sparskit matvec and external blas libs (true), don't use them (false)
  logical :: lilu  = .true.      ! use simple ilu preconditioner
 
  data epsmac/1.d-16/
 
  ! ilu0 preconditioner
 
  if (lilu) then
    ju0 = n+2
    jlu(1) = ju0 !!!
    iw = 0
    do ii = 1, n
      js = ju0
      do j=ia(ii),ia(ii+1)-1
        jcol = ja(j)
        if (jcol .eq. ii) then
          alu(ii) = aspar(j)
          iw(jcol) = ii
          ju(ii)  = ju0 !!!
        else
          alu(ju0) = aspar(j)
          jlu(ju0) = ja(j)
          iw(jcol) = ju0
          ju0 = ju0+1
        endif
      end do
      jlu(ii+1) = ju0 !!!
      jf = ju0-1
      jm = ju(ii)-1
      ! exit if diagonal element is reached.
      do j=js, jm
        jrow = jlu(j)
        tl = alu(j)*alu(jrow)
        alu(j) = tl
        ! perform  linear combination
        do jj = ju(jrow), jlu(jrow+1)-1
          jw = int(iw(jlu(jj)))
          if (jw .ne. 0) then
            alu(jw) = alu(jw) - tl*alu(jj)
            !                    write(*,*) ii, jw, jj
          end if
        end do
      end do
      !     invert  and store diagonal element.
      if (abs(alu(ii)) .lt. epsmac) then
        write (*,*) 'zero pivot'
        stop
      end if
      alu(ii) = 1.0d0/alu(ii)
      !     reset pointer iw to zero
      iw(ii) = 0
      do i = js, jf
        iw(jlu(i)) = 0
      end do
    end do
    ! end preconditioner
  end if
  !-------------------------------------------------------------
  its = 0
  ! outer loop starts here..
  if (lblas) then
    call amux (n, sol, vv, aspar, ja, ia)
  else
    do iii = 1, n
      t = 0.0d0
      do k = ia(iii), ia(iii+1)-1
        t = t + aspar(k) * sol(ja(k))
      end do
      vv(iii,1) = t
    end do
  end if
  do j=1,n
    vv(j,1) = rhs(j) - vv(j,1)
  end do
20 if (lblas) then
    ro = dnrm2(n, vv)
  else
    ro = sqrt(sum(vv(:,1)*vv(:,1)))
  end if
  if (abs(ro) .lt. epsmac) goto 999
  t = 1.0d0 / ro
  do j=1, n
    vv(j,1) = vv(j,1)*t
  end do
  if (its .eq. 0) eps1=eps*ro
  !      initialize 1-st term  of rhs of hessenberg system..
  rs(1) = ro
  i = 0
4 i=i+1
  its = its + 1
  i1 = i + 1
  if (lblas) then
    call lusol (n, vv(1,i), rhs, alu, jlu, ju)
    call amux (n, rhs, vv(1,i1), aspar, ja, ia)
  else
    do iii = 1, n !- lusol
      rhs(iii) = vv(iii,i)
      do k=jlu(iii),ju(iii)-1
        rhs(iii) = rhs(iii) - alu(k)* rhs(jlu(k))
      end do
    end do
    do iii = n, 1, -1
      do k=ju(iii),jlu(iii+1)-1
        rhs(iii) = rhs(iii) - alu(k)*rhs(jlu(k))
      end do
      rhs(iii) = alu(iii)*rhs(iii)
    end do
    do iii = 1, n !- amux
      t = 0.0d0
      do k = ia(iii), ia(iii+1)-1
        t = t + aspar(k) * rhs(ja(k))
      end do
      vv(iii,i1) = t
    end do
  end if
  !      modified gram - schmidt...
  if (lblas) then
    do j=1, i
      t = ddot(n, vv(1,j),vv(1,i1))
      hh(j,i) = t
      call daxpy(n, -t, vv(1,j), 1, vv(1,i1), 1)
      t = dnrm2(n, vv(1,i1))
    end do
  else
    do j=1, i
      t = 0.d0
      do iii = 1,n
        t = t + vv(iii,j)*vv(iii,i1)
      end do
      hh(j,i) = t
      vv(:,i1) = vv(:,i1) - t * vv(:,j)
      t = sqrt(sum(vv(:,i1)*vv(:,i1)))
    end do
  end if
  hh(i1,i) = t
  if ( abs(t) .lt. epsmac) goto 58
  t = 1.0d0/t
  do k=1,n
    vv(k,i1) = vv(k,i1)*t
  end do
  !     done with modified gram schimd and arnoldi step.. now  update factorization of hh
58 if (i .eq. 1) goto 121
  do k=2,i
    k1 = k-1
    t = hh(k1,i)
    hh(k1,i) = c(k1)*t + s(k1)*hh(k,i)
    hh(k,i) = -s(k1)*t + c(k1)*hh(k,i)
  end do
121 gam = sqrt(hh(i,i)**2 + hh(i1,i)**2)
  if (abs(gam) .lt. epsmac) gam = epsmac
  !      get  next plane rotation
  c(i) = hh(i,i)/gam
  s(i) = hh(i1,i)/gam
  rs(i1) = -s(i)*rs(i)
  rs(i) =  c(i)*rs(i)
  !      detrermine residual norm and test for convergence-
  hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i)
  ro = abs(rs(i1))
  if (i .lt. im .and. (ro .gt. eps1))  goto 4
  !      now compute solution. first solve upper triangular system.
  rs(i) = rs(i)/hh(i,i)
  do ii=2,i
    k=i-ii+1
    k1 = k+1
    t=rs(k)
    do j=k1,i
      t = t-hh(k,j)*rs(j)
    end do
    rs(k) = t/hh(k,k)
  end do
  !      form linear combination of v(*,i)'s to get solution
  t = rs(1)
  do k=1, n
    rhs(k) = vv(k,1)*t
  end do
  do j = 2, i
    t = rs(j)
    do k=1, n
      rhs(k) = rhs(k)+t*vv(k,j)
    end do
  end do
  !      call preconditioner.
  if (lblas) then
    call lusol (n, rhs, rhs, alu, jlu, ju)
  else
    do iii = 1, n
      do k=jlu(iii),ju(iii)-1
        rhs(iii) = rhs(iii) - alu(k)* rhs(jlu(k))
      end do
    end do
    do iii = n, 1, -1
      do k=ju(iii),jlu(iii+1)-1
        rhs(iii) = rhs(iii) - alu(k)*rhs(jlu(k))
      end do
      rhs(iii) = alu(iii)*rhs(iii)
    end do
  end if
  do k=1, n
    sol(k) = sol(k) + rhs(k)
  end do
  !      restart outer loop  when necessary
  if (ro .le. eps1) goto 990
  if (its .ge. maxits) goto 991
  !      else compute residual vector and continue..
  do j=1,i
    jj = i1-j+1
    rs(jj-1) = -s(jj-1)*rs(jj)
    rs(jj) = c(jj-1)*rs(jj)
  end do
  do j=1,i1
    t = rs(j)
    if (j .eq. 1)  t = t-1.0d0
    if (lblas) then
      call daxpy (n, t, vv(1,j), 1,  vv, 1)
    else
      vv(:,j) = vv(:,j) + t * vv(:,1)
    end if
  end do
  ! 199   format('   its =', i4, ' res. norm =', d20.6)
  !      restart outer loop.
  goto 20
990 ierr = 0
  return
991 ierr = 1
  return
999 continue
  ierr = -1
  return
  !---------------------------------------------------------------------

References amux(), daxpy(), ddot(), dnrm2(), and lusol().

qsplit()

subroutine qsplit	(	real*8, dimension(n)	a,
		integer, dimension(n)	ind,
		integer	n,
		integer	ncut
	)

Definition at line 3686 of file w3profsmd.F90.

  implicit none
 
  integer    :: n, ind(n), ncut
  real*8     :: a(n)
 
  !-----------------------------------------------------------------------
  !     does a quick-sort split of a real array.
  !     on input a(1:n). is a real array
  !     on output a(1:n) is permuted such that its elements satisfy:
  !
  !     abs(a(i)) .ge. abs(a(ncut)) for i .lt. ncut and
  !     abs(a(i)) .le. abs(a(ncut)) for i .gt. ncut
  !
  !     ind(1:n) is an integer array which permuted in the same way as a(*).
  !-----------------------------------------------------------------------
  real*8     :: tmp, abskey
  integer    :: itmp, first, last, j, mid
  !-----
  first = 1
  last = n
  if (ncut .lt. first .or. ncut .gt. last) return
  !
  !     outer loop -- while mid .ne. ncut do
  !
1 mid = first
  abskey = abs(a(mid))
  do j=first+1, last
    if (abs(a(j)) .gt. abskey) then
      mid = mid+1
      !     interchange
      tmp = a(mid)
      itmp = ind(mid)
      a(mid) = a(j)
      ind(mid) = ind(j)
      a(j)  = tmp
      ind(j) = itmp
    endif
  end do
  !
  !     interchange
  !
  tmp = a(mid)
  a(mid) = a(first)
  a(first)  = tmp
  !
  itmp = ind(mid)
  ind(mid) = ind(first)
  ind(first) = itmp
  !
  !     test for while loop
  !
  if (mid .eq. ncut) return
  if (mid .gt. ncut) then
    last = mid-1
  else
    first = mid+1
  endif
  goto 1
  !----------------end-of-qsplit------------------------------------------
  !-----------------------------------------------------------------------

Referenced by ilut().

runrc()

subroutine runrc	(	integer	n,
		real*8, dimension(n)	rhs,
		real*8, dimension(n)	sol,
		integer, dimension(16)	ipar,
		real*8, dimension(16)	fpar,
		real*8, dimension(*)	wk,
		real*8, dimension(n)	guess,
		real*8, dimension(*)	a,
		integer, dimension(*)	ja,
		integer, dimension(n+1)	ia,
		real*8, dimension(*)	au,
		integer, dimension(*)	jau,
		integer, dimension(*)	ju,
		external	solver
	)

Definition at line 3748 of file w3profsmd.F90.

  implicit none
  integer n,ipar(16),ia(n+1),ja(*),ju(*),jau(*)
  real*8 fpar(16),rhs(n),sol(n),guess(n),wk(*),a(*),au(*)
  external solver
  !-----------------------------------------------------------------------
  !     the actual tester. It starts the iterative linear system solvers
  !     with a initial guess suppied by the user.
  !
  !     The structure {au, jau, ju} is assumed to have the output from
  !     the ILU* routines in ilut.f.
  !
  !-----------------------------------------------------------------------
  !     local variables
  !
  integer       :: i, its
  !      real          :: dtime, dt(2), time
  !     external dtime
  save its
  !
  !     ipar(2) can be 0, 1, 2, please don't use 3
  !
  if (ipar(2).gt.2) then
    WRITE(*,*) 'I can not do both left and right preconditioning.'
    return
  endif
 
  its = 0
  !
  do i = 1, n
    sol(i) = guess(i)
  enddo
  !
  ipar(1) = 0
  !     time = dtime(dt)
 
10 call solver(n,rhs,sol,ipar,fpar,wk)
 
  if (ipar(7).ne.its) then
    its = ipar(7)
  endif
  if (ipar(1).eq.1) then
    call amux(n, wk(ipar(8)), wk(ipar(9)), a, ja, ia)
    goto 10
  else if (ipar(1).eq.2) then
    call atmux(n, wk(ipar(8)), wk(ipar(9)), a, ja, ia)
    goto 10
  else if (ipar(1).eq.3 .or. ipar(1).eq.5) then
    call lusol(n,wk(ipar(8)),wk(ipar(9)),au,jau,ju)
    goto 10
  else if (ipar(1).eq.4 .or. ipar(1).eq.6) then
    call lutsol(n,wk(ipar(8)),wk(ipar(9)),au,jau,ju)
    goto 10
  else if (ipar(1).le.0) then
    if (ipar(1).eq.0) then
      !            WRITE(*,*) 'Iterative sovler has satisfied convergence test.'
    else if (ipar(1).eq.-1) then
      WRITE(*,*) 'Iterative solver has iterated too many times.'
    else if (ipar(1).eq.-2) then
      WRITE(*,*) 'Iterative solver was not given enough work space.'
      WRITE(*,*) 'The work space should at least have ', ipar(4), &
           &           ' elements.'
    else if (ipar(1).eq.-3) then
      WRITE(*,*) 'Iterative sovler is facing a break-down.'
    else
      WRITE(*,*) 'Iterative solver terminated. code =', ipar(1)
    endif
  endif

References amux(), atmux(), lusol(), and lutsol().

Referenced by w3profsmd::w3xypfsnimp().

stopbis()

logical function stopbis	(	integer	n,
		integer, dimension(16)	ipar,
		integer	mvpi,
		real*8, dimension(16)	fpar,
		real*8, dimension(n)	r,
		real*8, dimension(n)	delx,
		real*8	sx
	)

Definition at line 2465 of file w3profsmd.F90.

  implicit none
  integer n,mvpi,ipar(16)
  real*8 fpar(16), r(n), delx(n), sx, ddot
  external ddot
  !-----------------------------------------------------------------------
  !     function for determining the stopping criteria. return value of
  !     true if the stopbis criteria is satisfied.
  !-----------------------------------------------------------------------
  if (ipar(11) .eq. 1) then
    stopbis = .true.
  else
    stopbis = .false.
  endif
  if (ipar(6).gt.0 .and. ipar(7).ge.ipar(6)) then
    ipar(1) = -1
    stopbis = .true.
  endif
  if (stopbis) return
  !
  !     computes errors
  !
  fpar(5) = sqrt(ddot(n,r,r))
  fpar(11) = fpar(11) + 2 * n
  if (ipar(3).lt.0) then
    !
    !     compute the change in the solution vector
    !
    fpar(6) = sx * sqrt(ddot(n,delx,delx))
    fpar(11) = fpar(11) + 2 * n
    if (ipar(7).lt.mvpi+mvpi+1) then
      !
      !     if this is the end of the first iteration, set fpar(3:4)
      !
      fpar(3) = fpar(6)
      if (ipar(3).eq.-1) then
        fpar(4) = fpar(1) * fpar(3) + fpar(2)
      endif
    endif
  else
    fpar(6) = fpar(5)
  endif
  !
  !     .. the test is struct this way so that when the value in fpar(6)
  !       is not a valid number, STOPBIS is set to .true.
  !
  if (fpar(6).gt.fpar(4)) then
    stopbis = .false.
    ipar(11) = 0
  else
    stopbis = .true.
    ipar(11) = 1
  endif
  !
  return

References ddot().

Referenced by bcgstab().

tidycg()

subroutine tidycg	(	integer	n,
		integer, dimension(16)	ipar,
		real*8, dimension(16)	fpar,
		real*8, dimension(n)	sol,
		real*8, dimension(n)	delx
	)

Definition at line 2523 of file w3profsmd.F90.

  implicit none
  integer i,n,ipar(16)
  real*8 fpar(16),sol(n),delx(n)
  !-----------------------------------------------------------------------
  !     Some common operations required before terminating the CG routines
  !-----------------------------------------------------------------------
  real*8 zero
  parameter(zero=0.0d0)
  !
  if (ipar(12).ne.0) then
    ipar(1) = ipar(12)
  else if (ipar(1).gt.0) then
    if ((ipar(3).eq.999 .and. ipar(11).eq.1) .or. &
         &        fpar(6).le.fpar(4)) then
      ipar(1) = 0
    else if (ipar(7).ge.ipar(6) .and. ipar(6).gt.0) then
      ipar(1) = -1
    else
      ipar(1) = -10
    endif
  endif
  if (fpar(3).gt.zero .and. fpar(6).gt.zero .and. ipar(7).gt.ipar(13)) then
    fpar(7) = log10(fpar(3) / fpar(6)) / dble(ipar(7)-ipar(13))
  else
    fpar(7) = zero
  endif
  do i = 1, n
    sol(i) = sol(i) + delx(i)
  enddo
  return

Referenced by bcgstab().

udsol()

subroutine udsol	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	au,
		integer, dimension(*)	jau
	)

Definition at line 3487 of file w3profsmd.F90.

  integer n, jau(*)
  real*8  x(n), y(n),au(*)
  !-----------------------------------------------------------------------
  !             Solves   U x = y  ;   U = upper triangular in MSR format
  !-----------------------------------------------------------------------
  ! solves a non-unit upper triangular matrix by standard (sequential )
  ! backward elimination - matrix stored in MSR format.
  ! with diagonal elements already inverted (otherwise do inversion,
  ! au(1:n) = 1.0/au(1:n),  before calling).
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real array containg the right side.
  !
  ! au,
  ! jau,    = Lower triangular matrix stored in modified sparse row
  !          format.
  !
  ! On return:
  !-----------
  !         x = The solution of  U x = y .
  !--------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8 t
  !-----------------------------------------------------------------------
  x(n) = y(n)*au(n)
  do k = n-1,1,-1
    t = y(k)
    do j = jau(k), jau(k+1)-1
      t = t - au(j)*x(jau(j))
    end do
    x(k) = au(k)*t
  end do
  !
  return
  !----------end-of-udsol-------------------------------------------------
  !-----------------------------------------------------------------------

udsolc()

subroutine udsolc	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	au,
		integer, dimension(*)	jau
	)

Definition at line 3575 of file w3profsmd.F90.

  integer n, jau(*)
  real*8 x(n), y(n), au(*)
  !-----------------------------------------------------------------------
  !    Solves     U x = y ;    U = nonunit Up. Triang. MSC format
  !-----------------------------------------------------------------------
  ! solves a (non-unit) upper triangular system by standard (sequential)
  ! forward elimination - matrix stored in Modified Sparse Column format
  ! with diagonal elements already inverted (otherwise do inversion,
  ! auuuul(1:n) = 1.0/au(1:n),  before calling ldsol).
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real*8 array containg the right hand side.
  !
  ! au,
  ! jau,   = Upper triangular matrix stored in Modified Sparse Column
  !          format.
  !
  ! On return:
  !-----------
  !         x = The solution of  U x = y .
  !--------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8 t
  !-----------------------------------------------------------------------
  do k=1,n
    x(k) = y(k)
  end do
  do  k = n,1,-1
    x(k) = x(k)*au(k)
    t = x(k)
    do j = jau(k), jau(k+1)-1
      x(jau(j)) = x(jau(j)) - t*au(j)
    end do
  end do
  !
  return
  !----------end-of-udsolc------------------------------------------------
  !-----------------------------------------------------------------------

uppdir()

subroutine uppdir	(	integer	n,
		real*8, dimension(n,lbp)	p,
		integer	np,
		integer	lbp,
		integer	indp,
		real*8, dimension(*)	y,
		real*8, dimension(*)	u,
		real*8, dimension(*)	usav,
		real*8	flops
	)

Definition at line 2384 of file w3profsmd.F90.

  implicit none
 
  integer   :: k,np,n,npm1,j,ju,indp,lbp
  real*8    :: p(n,lbp), y(*), u(*), usav(*), x, flops
 
  !-----------------------------------------------------------------------
  !     updates the conjugate directions p given the upper part of the
  !     banded upper triangular matrix u.  u contains the non zero
  !     elements of the column of the triangular matrix..
  !-----------------------------------------------------------------------
  real*8 zero
  parameter(zero=0.0d0)
  !
  npm1=np-1
  if (np .le. 1) goto 12
  j=indp
  ju = npm1
10 if (j .le. 0) j=lbp
  x = u(ju) /usav(j)
  if (x .eq. zero) goto 115
  do k=1,n
    y(k) = y(k) - x*p(k,j)
  end do
  flops = flops + 2*n
115 j = j-1
  ju = ju -1
  if (ju .ge. 1) goto 10
12 indp = indp + 1
  if (indp .gt. lbp) indp = 1
  usav(indp) = u(np)
  do  k=1,n
    p(k,indp) = y(k)
  end do
  return
  !-----------------------------------------------------------------------
  !-------end-of-uppdir---------------------------------------------------

usol()

subroutine usol	(	integer	n,
		real*8, dimension(n)	x,
		real*8, dimension(n)	y,
		real*8, dimension(*)	au,
		integer, dimension(*)	jau,
		integer, dimension(n+1)	iau
	)

Definition at line 3444 of file w3profsmd.F90.

  integer n, jau(*),iau(n+1)
  real*8  x(n), y(n), au(*)
  !-----------------------------------------------------------------------
  !             Solves   U x = y    U = unit upper triangular.
  !-----------------------------------------------------------------------
  ! solves a unit upper triangular system by standard (sequential )
  ! backward elimination - matrix stored in CSR format.
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real array containg the right side.
  !
  ! au,
  ! jau,
  ! iau,    = Lower triangular matrix stored in compressed sparse row
  !          format.
  !
  ! On return:
  !-----------
  !         x = The solution of  U x = y .
  !--------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8  t
  !-----------------------------------------------------------------------
  x(n) = y(n)
  do  k = n-1,1,-1
    t = y(k)
    do j = iau(k), iau(k+1)-1
      t = t - au(j)*x(jau(j))
    end do
    x(k) = t
  end do
  !
  return
  !----------end-of-usol--------------------------------------------------
  !-----------------------------------------------------------------------

usolc()

subroutine usolc	(	integer	n,
		real*8, dimension(*)	x,
		real*8, dimension(*)	y,
		real*8, dimension(*)	au,
		integer, dimension(*)	jau,
		integer, dimension(*)	iau
	)

Definition at line 3531 of file w3profsmd.F90.

  real*8  x(*), y(*), au(*)
  integer n, jau(*),iau(*)
  !-----------------------------------------------------------------------
  !       SOUVES     U x = y ;    where U = unit upper trang. CSC format
  !-----------------------------------------------------------------------
  ! solves a unit upper triangular system by standard (sequential )
  ! forward elimination - matrix stored in CSC format.
  !-----------------------------------------------------------------------
  !
  ! On entry:
  !----------
  ! n      = integer. dimension of problem.
  ! y      = real*8 array containg the right side.
  !
  ! au,
  ! jau,
  ! iau,    = Uower triangular matrix stored in compressed sparse column
  !          format.
  !
  ! On return:
  !-----------
  !         x  = The solution of  U x  = y.
  !-----------------------------------------------------------------------
  ! local variables
  !
  integer k, j
  real*8 t
  !-----------------------------------------------------------------------
  do k=1,n
    x(k) = y(k)
  end do
  do  k = n,1,-1
    t = x(k)
    do j = iau(k), iau(k+1)-1
      x(jau(j)) = x(jau(j)) - t*au(j)
    end do
  end do
  !
  return
  !----------end-of-usolc-------------------------------------------------
  !-----------------------------------------------------------------------

vbrmv()

subroutine vbrmv	(	integer	nr,
		integer	nc,
		integer, dimension(nr+1)	ia,
		integer, dimension(*)	ja,
		integer, dimension(*)	ka,
		real*8, dimension(*)	a,
		integer, dimension(nr+1)	kvstr,
		integer, dimension(*)	kvstc,
		real*8, dimension(*)	x,
		real*8, dimension(*)	b
	)

Definition at line 3161 of file w3profsmd.F90.

  !-----------------------------------------------------------------------
  integer nr, nc, ia(nr+1), ja(*), ka(*), kvstr(nr+1), kvstc(*)
  real*8  a(*), x(*), b(*)
  !-----------------------------------------------------------------------
  !     Sparse matrix-full vector product, in VBR format.
  !-----------------------------------------------------------------------
  !     On entry:
  !--------------
  !     nr, nc  = number of block rows and columns in matrix A
  !     ia,ja,ka,a,kvstr,kvstc = matrix A in variable block row format
  !     x       = multiplier vector in full format
  !
  !     On return:
  !---------------
  !     b = product of matrix A times vector x in full format
  !
  !     Algorithm:
  !---------------
  !     Perform multiplication by traversing a in order.
  !
  !-----------------------------------------------------------------------
  !-----local variables
  integer n, i, j, ii, jj, k, istart, istop
  real*8  xjj
  !---------------------------------
  n = kvstc(nc+1)-1
  do i = 1, n
    b(i) = 0.d0
  enddo
  !---------------------------------
  k = 1
  do i = 1, nr
    istart = kvstr(i)
    istop  = kvstr(i+1)-1
    do j = ia(i), ia(i+1)-1
      do jj = kvstc(ja(j)), kvstc(ja(j)+1)-1
        xjj = x(jj)
        do ii = istart, istop
          b(ii) = b(ii) + xjj*a(k)
          k = k + 1
        enddo
      enddo
    enddo
  enddo
  !---------------------------------
  return