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

Popular posts from this blog

How to provide Authorization & Authentication using Asp.net, C#? -

toolbar - How to add link to user registration inside toobar in admin joomla 3 custom component -

How to use Authorization & Authentication in Asp.net, C#? -