Calling Numerical Recipes routines can be tricky. We need to make sure that the parameters are passed correctly to the subroutines/functions. Numerical Recipes makes the matters worse because of the poor programming. For example, in the GAUSSJ routine, the parameter mp is not the leading dimension of B -- it is the second dimension of B -- and is not needed. Because Numerical Recipes uses real numbers as matrix/vector indices (which is a sloppy programming practice as index variables should be integers), be sure to declare dimensions with real numbers. Actually Numerical Recipes is not the only one to blame. Bill Gates and his QBasic does not make clear distinction of real numbers and integers, which comes back to haunt you poor students. Now is the payback time for being sloppy. Below we show some examples that work and some examples that are no good. To make matters even worse, the buggy program runs and yields results, but they are wrong and highly unpredictable. They happen without any warning -- no complaints of syntax errors from compilers or linkers nor any indication during run-time!!! And you may not be aware that the answers are wrong unless you test them thoroughly.
DECLARE SUB gaussj(a!(), n!, np!, b!(), m!, mp!)
const n=3 ' n is treated as an integer
dim A(n,n), b(n,1)
:
call GAUSSJ(A(), n, n, b(), 1, 1)
The rest of the codes for an example of A*x=b are in:
DECLARE SUB gaussj(a!(), n!, np!, b!(), m!, mp!)
const n=3! ' n is treated as a real number
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1, 1)
Likewise, the following is o.k., because the period "." forces n
to be a real number
DECLARE SUB gaussj(a!(), n!, np!, b!(), m!, mp!)
const n=3. ' n is treated as a real number
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1, 1)
The following is no good, because the leading dimension of b
(i.e., the constant "n") is messed up, in spite of the
declaration of SUB and what type of arguments the subroutine
should take. Moral: DECLARE alone cannot ensure that the
variable types match, even though DECLARE was meant to catch
type mismatch.
DECLARE SUB gaussj(a!(), n!, np!, b!(), m!, mp!)
const n=3 ' n is treated as an integer
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1, 1)
Likewise, the following is no good, because it has the same
problem with "n" as the last example, although the "1!" part has been
fixed.
DECLARE SUB gaussj(a!(), n!, np!, b!(), m!, mp!)
const n=3 ' n is treated as an integer
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1!, 1!)
The one below is o.k. because the problem with "n" is fixed.
DECLARE SUB gaussj(a!(), n!, np!, b!(), m!, mp!)
n=3 ' n is treated as a real variable by default
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1!, 1!)
The following is o.k., because all the arguments are real
const n=3! ' n is treated as a real number
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1!, 1!)
The following is also o.k.
n=3 ' n is treated as a real variable by default
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1!, 1!)
The following is also o.k.
n=3. ' n is treated as a real variable by default
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1., 1.)
The following is no good because n is treated as an integer
const n=3 ' n is treated as an integer
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1!, 1!)
The following is no good because 1 is treated as an integer
n=3 ' n is treated as a real variable by default
dim A(n,n), b(n)
:
call GAUSSJ(A(), n, n, b(), 1, 1)
There are more permutations that are not listed above, but I think we have had enough. It is easy to get confused after all these different permutations. The bottom line is that you want to distinguish clearly what is an integer and what is a real number.
|
