[ Pobierz całość w formacie PDF ]

partitioned into blocks of sizes as shown below:
îø ùø
A : m × n RHS : m × p
ïø úø
ïø úø
C = . (3.3.1)
ðø ûø
Ä : 1 × n 0 : 1 × p
Now a good way to begin the writing of a program such as the general-purpose matrix
94 Numerical linear algebra
analysis program that we now have in mind is to consider the different procedures, or
modules, into which it may be broken up. We suggest that the individual blocks that we
are about to discuss should be written as separate subroutines, each with its own clearly
defined input and output, each with its own documentation, and each with its own local
variable names. They should then be tested one at a time, by giving them small, suitable
test problems. If this is done, then the main routine won t be much more than a string of
calls to the various blocks.
1. Proceduresearchmat(C,r,s,i1,j1,i2,j2)
This routine is given an r × s array C, and two positions in the matrix, say [i1, j1] and
[i2, j2]. It then carries out a search of the rectangular submatrix of C whose northwest
corner is at [i1, j1] and whose southeast corner is at [i2, j2], inclusive, in order to find an
element of largest absolute value that lives in that rectangle. The subroutine returns this
element of largest magnitude, asbig, and the row and column in which it lives, asiloc,
jloc.
Subroutinesearchmat will be called in at least two different places in the main routine.
First, it will do the search for the next pivot element in the southeast rectangle. Second, it
can be used to determine if the equations are consistent by searching the right-hand sides
of equations r +1, . . . , m (r is the pseudorank) to see if they are all zero (i.e., below our
tolerance level).
2. Procedureswitchrow(C,r,s,i,j,k,l)
The program is given an r × s matrix C, and four integers i, j, k and l. The subroutine
interchanges rows i and j of C, between columns k and l inclusive, and returns a variable
calledsign with a value of -1, unless i = j, in which case it does nothing to C and returns
a +1 insign.
3. Procedureswitchcol(C,r,s,i,j,k,l)
This subroutine is like the previous one, except it interchanges columns i and j of C,
between rows k and l inclusive. It also returns a variable calledsign with a value of -1,
unless i = j, in which case it does nothing to C and returns a +1 insign.
The subroutinesswitchrow andswitchcol are used during the forward solution in the
obvious way, and again after the back solution has been done, to unscramble the output
(see procedure 5 below).
4. Procedurepivot(C,r,s,i,k,u)
Given the r × s matrix C, and three integers i, k and u, the subroutine assumes that
Cii =1. It then stores Cki in the local variable tm sets Cki to zero, and reduces row k of
C, in columns u to s, by doing the operation Ckq := Ckq - t " Ciq for q = u, . . . , s.
The use of the parameter u in this subroutine allows the flexibility for economical op-
eration in both the forward and back solution. In the forward solution, we take u = i +1
and it reduces the whole row k. In the back solution we use u = n + 1, because the rest of
row k will have already been reduced.
5. Procedurescale(C,r,s,i,u)
3.3 Building blocks for the linear equation solver 95
Given an r × s matrix C and integers i and u, the routine stores Cii in a variable called
piv. It then does Cij := Cij/piv for j = u, . . . , s and returns the value of piv.
6. Procedurescramb(C,r,s,n)
This procedure permutes the first n rows of the r × s matrix C according to the permu-
tation that occupies the positions Cr1, Cr2, . . . , Crn on input.
The use of this subroutine is explained in detail in section 3.6 (q.v.). Its purpose is to
rearrange the rows of the output matrix that holds a basis for the kernel, and also the rows of
the output matrix that holds particular solutions of the give system(s). After rearrangement
the rows will correspond to the original numbering of the unknowns, thereby compensating
for the renumbering that was induced by column interchanges during the forward solution.
This subroutine poses some interesting questions if we require that it should not use any
additional array space beyond the input matrix itself.
7. Procedureident(C,r,s,i,j,n,q)
This procedure inserts q times the n × n identity matrix into the n × n submatrix whose
Northwest corner is at position [i, j] of the r × s matrix C.
Now let s look at the assembly of these building blocks into a complete matrix analysis
procedure calledmatalg(C,r,s,m,n,p,opt,eps). Input items to it are:
" An r ×s matrix C (as well as the values of r and s), whose Northwest m×n submatrix
contains the matrix of coefficients of the system(s) of equations that are about to be
solved. The values of m and n must also be provided to the procedure. It is assumed
that r =1 + max(m, n). Unless the inverse of the coefficient matrix is wanted, the
Northeast m × p submatrix of C holds p different right-hand side vectors for which
we want solutions.
" The numbers r, s, m, n and p.
" A parameteroption that is equal to 1 if we want an inverse, equal to 2 if we want
to see the determinant of the coefficient matrix (if square) as well as a basis for the
kernel (if it is nontrivial) and a set of p particular solution vectors.
" A real parametereps that is used to bound roundoff error.
Output items from the procedurematalg are:
" The pseudorank r
" The determinantdet if m = n
" An n × r matrixbasis , whose columns are a basis for the kernel of the coefficient
matrix.
" An n × p matrixpartic, whose columns are particular solution vectors for the given
systems.
96 Numerical linear algebra
In caseopt = 1 is chosen, the procedure will fill the last m columns and rows of C with
an m × m identity matrix, set p = n = m, and proceed as before, leaving the inverse matrix
in the same place, on output.
Let s remark on how the determinant is calculated. The reduction of the input matrix [ Pobierz całość w formacie PDF ]