How are these arrays being used in this Fortran algorithm? -
i'm writing subroutines in fortran90 perform numerical computations. however, part of this, need include codes netlib templates library written in fortran77. i'm having hard time getting them work - understanding usage of arrays.
for instance, need use subroutine called gmres. here arguments:
subroutine gmres( n, b, x, restrt, work, ldw, work2, ldw2, $ iter, resid, matvec, psolve, info ) * .. scalar arguments .. integer n, restrt, ldw, ldw2, iter, info double precision resid * .. * .. array arguments .. double precision b( * ), x( * ), work( * ), work2( * ) * .. * .. function arguments .. * external matvec, psolve * * * arguments * ========= * * n (input) integer. * on entry, dimension of matrix. * unchanged on exit. * * b (input) double precision array, dimension n. * on entry, right hand side vector b. * unchanged on exit. * * x (input/output) double precision array, dimension n. * on input, initial guess; on exit, iterated solution. * * restrt (input) integer * restart parameter, <= n. parameter controls amount * of memory required matrix work2. * * work (workspace) double precision array, dimension (ldw,5). * note if initial guess 0 vector, * storing initial residual not necessary. * * ldw (input) integer * leading dimension of array work. ldw >= max(1,n). * * work2 (workspace) double precision array, dimension (ldw2,2*restrt+2). * workspace used constructing , storing * upper hessenberg matrix. 2 columns used * store givens rotation matrices. * * ldw2 (input) integer * leading dimension of array work2. * ldw2 >= max(1,restrt). * * iter (input/output) integer * on input, maximum iterations performed. * on output, actual number of iterations performed. * * resid (input/output) double precision * on input, allowable error tolerance. * on ouput, norm of residual vector if solution * approximated tolerance, otherwise reset input * tolerance. * * info (output) integer * = 0: successful exit * = 1: maximum number of iterations performed; * convergence not achieved. * * matvec (external subroutine) * user must provide subroutine perform * matrix-vector product a*x = y. * vector x must remain unchanged. solution * over-written on vector y. * * call is: * * call matvec( x, y ) notice how arrays defined work( * ), work2( * ). mind these 1d arrays of arbitrary length. in argument description, seems suggesting 2d arrays, matrices, dimension work(ldw, 5). 1d or 2d?
also, in gmres algorithm, these work arrays used this:
call matvec(sclr1, work(ndx1), sclr2, work(ndx2)) so if work arrays 2d, wouldn't above give rank mismatch? mean access 2d array 1 index that? should define work arrays 2d or 1d?
edit
the gmres routine requires matvec routine supplied. in gmres code being called
call matvec(sclr1, x, sclr2, work(ndx2)) and like
call matvec(sclr1, work(ndx1), sclr2, work(ndx2)) my subroutine matvec i'm trying supply looks like:
subroutine matvec(alpha, x, beta, y) real(dp), intent(in) :: alpha, beta real(dp), dimension(*), intent(in) :: x real(dp), dimension(*), intent(inout) :: y real(dp), dimension(*,*) :: integer :: n call make_jac(n,a) call dgemv('notranspose', n, n, alpha, a, n, x, 1, beta, y, 1) end subroutine matvec where make_jac returns matrix problem i'm working on, along dimension n. blas routine dgemv handles matrix-vector product.
work( * ) declares assumed size array, can two-dimensional. see here.
the compiler not complain if feed one-dimensional array subroutine, weird things (up , including segmentation fault) might happen.
better use arrays matching specifications.
Comments
Post a Comment