From the top level of the Help navigator, you can find a summary of IDL's linear systems routines by using the "Navigate" button followed by "Chapter Topics" then "List of Routines by Application" then "Mathematics Routines" and then "Linear Systems Routines". Alternately, you can find a more detailed description using the "Navigate" button and selecting "Contents" and then the Help book "Using IDL", chapter "Mathematics", section "Linear Systems".
IDL is able to manipulate arrays (for example, matrices and vectors) by referring to the array descriptor. There is also a concise notation for matrix multiplication: a##b multiplies two matrices (a and b) together (see Language Features Quick Reference/Operators).
Here is how you set up two matrices and multiply them together using IDL:
IDL> array1 = [ [1, 2, 1], [2, -1, 2] ]
IDL> print,array1
1 2 1
2 -1 2
IDL> array2 = [ [1, 3], [0, 1], [1, 1] ]
IDL> print,array2
1 3
0 1
1 1
IDL> PRINT, array1##array2
2 6
4 7
All this seems to be fine. But now try
the following matrix/vector multiplies:
IDL> array3=[1, 0, 1] ; This should appear as a row vector
IDL> print,array3 ; Yup. This is a row alright.
1 0 1
IDL> print,array1##array3 ;OOPS! array3 is considered to be a column!
2
4
IDL> array4=[ [1], [0], [1] ] ; This should appear as a column vector!
IDL> print,array4 ; Yup. This is a column.
1
0
1
IDL> print,array1##array4 ; Hmmm. So array3 = array 4?
2
4
IDL> print,array3-array4
0 0 0
IDL> print,array4-array3
0
0
0
IDL> print,array4##array3 ; Now this will test it.......!
1 0 1
0 0 0
1 0 1
IDL> print,array3##array4
2
So vectors (1-D arrays) appear to be interpreted
as column vectors or row vectors depending on the
context.
Now look at the following incompatible system:
IDL> array5=[[1, 0, 1], [3, 1, 1]]
IDL> print,array5
1 0 1
3 1 1
IDL> print,array1##array5
% Operands of matrix multiply have incompatible dimensions: ARRAY5, ARRAY1.
% Execution halted at: $MAIN$
Once again, (2D) arrays are interpreted as they look!
The routine 'ludc' finds the LU decomposition of a square matrix with partial pivoting. You follow it with a call to 'lusol' to perform the back substitution and solve the system. Both of these routines (and many others from the Linear Systems menu) are taken from the book "Numerical Recipes". Let's follow the following example:
IDL> a = [[ 2.0, 1.0, 1.0], [ 4.0, -6.0, 0.0], [-2.0, 7.0, 2.0]]
IDL> print,a ;coefficient matrix
2.00000 1.00000 1.00000
4.00000 -6.00000 0.00000
-2.00000 7.00000 2.00000
IDL> atemp=a ;save this for later (as a is overwritten)
IDL> b= [[3.0], [-8.0], [10.0]] & print,b ;column right hand side??
3.00000
-8.00000
10.0000
IDL> ludc,a,index ;LU decomposition of a
IDL> result = lusol(a,index,b) ;back substitution
% LUSOL: Argument B must be a column vector: B
% Execution halted at: $MAIN$
IDL> b = [3.0, -8.0, 10.] & print,b ; NO. The documentation says that
b must be a VECTOR!!
3.00000 -8.00000 10.0000
IDL> result = lusol(a,index,b)
IDL> print,result ;This works!
1.00000 2.00000 -1.00000
IDL> print,atemp##result ;R.H.S. is printed as a column vector!
3.00000
-8.00000
10.0000
So all this seems to do the job. Just remember that a "vector" is a 1D array,
not an array of 1-element arrays.
I would suggest that whenever you
wish to solve a big matrix equation in IDL that you experiment with a
small system first to make sure that you have the conventions right.
IDL> a=findgen(4,4)*0. + 1.&print,a
1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000 1.00000
1.00000 1.00000 1.00000 1.00000
creates a unit floating point array. You can also call the
"make_array" subroutine.