Let's look at a specific example. Let's assume you want to use ELL storage with with a maximum of 7 non-zeros per row, that the largest systems you want to solve have 1000 unknowns, and that this particular system has 800 unknowns. In this case you'd have the following (or something similar) in your main program:
integer ia(_JT_no_of_storage_parameters_) ! storage info for coefficient
integer iap(_JT_no_of_storage_parameters_) ! storage info for preconditioner
integer ja(1000,7) ! column map array for coefficient
integer jap(1000,7) ! column map array for preconditioner
real a(1000,7) ! coefficient matrix
real ap(1000,7) ! preconditioner
real x(1000), b(1000) ! solution and source vectors
Note that I've used a sans serif font for the cpp variables to emphasize the fact that they will be replaced by their appropriate definitions by cpp.
Then set up array ia as follows:
ia(_JT_storage_) = _JT_storage_ELL_
ia(_JT_nunk_) = 800
ia(_JT_leading_dimension_) = 1000
ia(_JT_maxnz_) = 7
iap(_JT_storage_) = _JT_storage_ELL_
iap(_JT_nunk_) = 800
iap(_JT_leading_dimension_) = 1000
iap(_JT_maxnz_) = 7
I've omitted comments, since use of cpp variables makes the code virtually self-documenting.
The remaining elements of ia would not be used, and need not be set.
Note that ja is the same shape as a (has the same dimensions as a). For ELL storage, ja basically contains the true column position of each element in a (see the example above in the discussion of the elements of ia). If a different storage that did not require ja, such as full conventional storage, were used, ja could simply be a scalar.
In addition, you need vectors for the integer and real control parameters:
integer iparm(_JT_no_of_iparms_)
real rparm(_JT_no_of_rparms_)
cpu times:
real cpu(0:500)
for norms of the estimates of the residual and error:
real rnorm(0:500)
real err(0:500)
and for the norms of the true residual and error:
real rnormt(0:500)
real errt(0:500)
Here I've set the upper dimension of the residual and error norm arrays to 500. Note that iparm(_JT_ itmax_) must be set less than or equal to this value.
Since rnormt and errt are only calculated if desired, if they are not calculated they can be scalars. However, bad things will certainly happen if you make them scalars and then later decide you'd like to calculate them and forget to give them adequate storage space.
Next you would call JT_ SETDEFAULTS in your main program to set default values for iparm and rparm:
call JT_SetDefaults (iparm, rparm, status)
This sets reasonable values for several of the parameters, such as which norm to use, the maximum number of iterations to allow, convergence criterion, etc. The default values are noted above. Then you can set particular values manually as you wish.
After setting up the coefficient and source, solve the system by calling JT_ GMRES as follows:
call JT_GMRES (b, a, ia, ja, iparm, rparm, x,
& cpu, rnormt, errt, rnorm, err, ap, iap, jap, status)
which should be followed by code to check the return status. For example:
if (status .eq. -1) then
write(iparm(4),*) 'invalid arguments to JT_GMRES'
[stop, return, etc.]
elseif (status .eq. -2) then
write(iparm(4),*) 'memory allocation failure within JT_GMRES'
[stop, return, try a smaller problem, etc.]
elseif (status .eq. -3) then
write(iparm(4),*) 'internal error in JT_GMRES'
[stop, return, etc.]
elseif (status .eq. -4) then
write(iparm(4),*) 'itmax exceeded in JT_GMRES'
[stop, return, restart, try a different solver, etc.]
elseif (status .eq. -5) then
write(iparm(4),*) 'breakdown in iterates in JT_GMRES'
[stop, return, try a different solver, etc.]
endif
One thing to note is that if itmax is exceeded or if a breakdown in iterates occurs, one should probably examine the residual norm and error estimate and consider accepting the solution obtained, since it might be ``good enough''.