Each process(or) p gets lmatbk(p) consecutive rows of the GLOBAL matrix A starting at row ptrmbk(p) +1. Thus we have a LOCAL lmatbk(p) x l matrix A(p) on every process(or) p . This LOCAL matrix is stored in the sparse matrix storage scheme described in the next chapter. Additionaly, LINSOL assumes a logical column block splitting of the matrix on each process(or). This is explained in the chapter Column Block Splitting .
A(p) = A(1,p) + A(2,p) + ... + A(i,p) + ... + A(10,p)
Each submatrix A(i,p) contains one (and only one) special data structure. The vectors within the allowed data structures are named "vector terms". The number of vector terms for A(i,p) is stored in the (virtual) array nvt(i,p) . The different submatrices (and vector terms respectively) can be assembled without cutback to the matrix A(p) .
Two elements of two different vector terms can have the same entry in the matrix A(p) - this means that the entry is the sum of the two elements. The overlay of elements in the matrix A(p) does not influence the correct computation of the solution of the linear system; but then the normalization does not work in the specified way - this means that you don't get exactly the chosen normalization (except of the normalization by main diagonal).
In the following, the index p indicating the process(or) number will be dropped. The integer number typ indicates the type of the LOCAL submatrix.
The variables iac , iar and the array index are GLOBAL and relate to the GLOBAL matrix A . The variables l and lvt are local and relate to the LOCAL matrix A . The variables adda , indc and indr are pointers and point to the matrix array mat and to the array of indices index respectively. The following examples show the settings of the varables for the first block of the LOCAL matrix on the first process(or)- this means that for the other blocks and (or) on the other processes (processors) the values of the GLOBAL variables ( iac , iar , index ) change.
MAIN : main diagonal ( typ =10 or 1)
A(1) is the main diagonal of the matrix A . This is the diagonal starting with row- and column-index ptrmbk(p) +1. Gaps with no entries have to be filled by zeroes
1234567890123... 1 *------- 2 -*------ 3 --*----- A(1) = 4 ---*---- 5 ----*--- 6 -----*-- 7 ------*- 8 -------*
FULLDIAG : full diagonal ( typ =20 or 2)
Starting at the column iac+1 and the row iar+1 , the vector term FULLDIAG is a diagonal portion of length lvt . Gaps with no entries have to be filled by zeroes.
In the example iac =3, iar =1 and lvt =4 holds.
1234567890123... 1 -------- 2 ---*---- 3 ----*--- A(2) = 4 -----*-- 5 ------*- 6 -------- 7 -------- 8 --------
PDIAG : packed diagonal ( typ =30 or 3)
Starting at the column iac+1 and the row iar+1 , the vector term PDIAG is a sparse diagonal portion with length lvt . A pointer indc points to the starting adress of the array index which is used to adress the nonzero elements in the diagonal.
In the example iac =3, iar =1, lvt =3 and index(indc+1..lvt) =(1,2,4) holds.
1234567890123... 1 -------- 2 ---*---- 3 ----*--- A(3) = 4 -------- 5 ------*- 6 -------- 7 -------- 8 --------
IXCOL : indexed column ( typ =40 or 4)
This vector term has entries in lvt contiguous rows starting at row iar+1 . A pointer indc points to the starting adress of the array index which is used to adress the column indices relating to the contiguous row indices.
In the example iar =1, lvt =4 and index(indc+1..lvt) =(4,2,1,7) holds.
1234567890123... 1 -------- 2 ---*---- 3 -*------ A(4) = 4 *------- 5 ------*- 6 -------- 7 -------- 8 --------
IXROW : indexed row ( typ =50 or 5)
This vector term has entries in lvt contiguous columns starting at column iac+1 . A pointer indr points to the starting adress of the array index which is used to adress the row indices relating to the contiguous column indices.
In the example iac =0, lvt =4 and index(indr+1..lvt) =(4,3,6,1) holds.
1234567890123... 1 ---*---- 2 -------- 3 -*------ A(5) = 4 *------- 5 -------- 6 --*----- 7 -------- 8 --------
SKY : starry sky ( typ =60 or 6)
lvt entries of nonzero elements are completely defined by two pointers indc,indr . The column-indices are listed in the integer vector index(indc+1..lvt) , the row-indices in the vector index(indr+1..lvt) .
In the example lvt =3, index(indc+1..lvt) =(1,2,6), and index(indr+1..lvt) =(3,1,7) holds.
1234567890123... 1 -*------ 2 -------- 3 *------- A(6) = 4 -------- 5 -------- 6 -------- 7 -----*-- 8 --------
FULLROW: full row ( typ =70 or 7)
FULLROW defines lvt contiguous elements in the row iar +1 starting at column iac +1. Gaps with no entries have to be filled by zeroes.
In the example iar =7, iac =2 and lvt =5 holds.
1234567890123... 1 -------- 2 -------- 3 -------- A(7) = 4 -------- 5 -------- 6 -------- 7 -------- 8 --*****-
FULLCOL : full column ( typ =80 or 8)
FULLCOL defines lvt contiguous elements in the column iac +1 starting at row iar +1. Gaps with no entries have to be filled by zeroes.
In the example iac =7, iar =2 and lvt =3 holds.
1234567890123... 1 -------- 2 -------- 3 -------* A(8) = 4 -------* 5 -------* 6 -------- 7 -------- 8 --------
PROW: packed row ( typ =90 or 9) {not yet implemented!!!!}
PROW defines lvt elements adressed by vector index via indc in the row iar +1.
In the example iar =7, lvt =3 and index(indc+1..lvt) =(3,5,8) holds.
1234567890123... 1 -------- 2 -------- 3 -------- A(9) = 4 -------- 5 -------- 6 -------- 7 -------- 8 --*-*--*
PCOL : packed column ( typ =100 or 11) {not yet implemented!!!!}
PCOL defines lvt elements adressed by vector index via indr in the column iac +1.
In the example iac =7, lvt =3 and index(indr+1..lvt) =(3,4,7) holds.
1234567890123... 1 -------- 2 -------- 3 -------* A(11) = 4 -------* 5 -------- 6 -------- 7 -------* 8 --------
On each process, all index vectors are stored in the integer array
index
of length
lindex
.
Indc
+1 and
indr
+1 point to the first
element of index vectors relating to vector terms within the array
index
.
The entries of the LOCAL matrix
A
are stored in the real array
mat
of length
lmat
.
Adda
+1 points to the first element
of a vector term within the array
mat
.
The integer numbers
typ, adda, lvt, iac, iar, indc
and
indr
are
handed over in the integer matrix
info(ia1,ia2)
;
ia1>=nvt
and
ia2
>=7 must hold. The values
info(i,1)
,...,
info(i,7)
specify the information for the i-th vector term of
typ=info(i,1)
.
The following table shows the meaning of
info(i,k)
for the different data structures:
shape | k=1 | k=2 | k=3 | k=4 | k=5 | k=6 | k=7 |
MAIN | 1 or 10 | adda | l | 0 | 0 | 0 | 0 |
FULLDIAG | 2 or 20 | adda | lvt | iac | iar | 0 | 0 |
PDIAG | 3 or 30 | adda | lvt | iac | iar | indc | 0 |
IXCOL | 4 or 40 | adda | lvt | 0 | iar | indc | 0 |
IXROW | 5 or 50 | adda | lvt | iac | 0 | 0 | indr |
SKY | 6 or 60 | adda | lvt | 0 | 0 | indc | indr |
FULLROW | 7 or 70 | adda | lvt | iac | iar | 0 | 0 |
FULLCOL | 8 or 80 | adda | lvt | iac | iar | 0 | 0 |
PROW | 9 or 90 | adda | lvt | 0 | iar | indc | 0 |
PCOL | 11 or 100 | adda | lvt | iac | 0 | 0 indr |
A = AA + D + AA(T)
where D is the main diagonal of A and AA is a M X M -matrix (not necessarily equal to the upper or lower triangle of A ). AA(T) denotes the transposed matrix of AA . In this case only the matrix AA and the main diagonal (if nonzero) must be handed over.First we look to a nonsymmetric problem. MAT is a 8 X 8 matrix (=>l=8):
| 4. 0 2. 0 0 0 5. 0. | | 0 4. -1. 3. 0 0 0 0. | | 0 0 4. -1. 0 -2. 0 0. | MAT = | 3. 0 0 4. -1. 0 0 0. | | 0 0 0 0 4. -1. 3. 0. | | 0 0 0 0 0 4. 0 0. | | 0 0 0 0 0 0 4. 0. | | 0 0 1. 1. 1. 1. 1. 4. |
The matrix MAT is subdivided into five vector terms
MAT = MAIN + 1 x FULLDIAG + 1 x IXCOL + 1 x SKY + 1 x FULLROW = A(1) + A(2) + A(4) + A(6) + A(7) .We set ia1=nvt=5 and ia2=7. We get
| 4. | | 4. | | 4. | A(1) = | 4. | | 4. | | 4. | | 4. | | 4.|
| | | -1. | | -1. | A(2) = | -1. | | -1. | | | | | | |
| 5. | | 3. | | -2. | A(4) = | 3. | | | | | | | | |
| 2. | | | | | A(6) = | | | 3. | | | | | | |
| | | | | | A(7) = | | | | | | | | | 1. 1. 1. 1. 1. | .
[------A(1)-----------][----A(2)------][----A(4)---][-A(6)][-----A(7)-----] MAT = (4.,4.,4.,4.,4.,4.,4.,4.,-1.,-1.,-1.,-1.,5.,3.,-2.,3.,2.,3.,1.,1.,1.,1.,1.) ^ ^ ^ ^ ^ adda(1)=0 adda(2)=8 adda(3)=12 adda(4)=16 adda(5)=18The index vectors are stored in the array index (lindex=8):
[--A(4)-][--A(6)--] index = (7,4,6,1,3,7, 1,5) ^ ^ ^ info(3,6)=0 info(4,4)=4 info(4,6)=6
| 4. 0 2. 3. 0 0 5. 0. | | 0 4. -1. 3. 0 0 0 0. | | 2 -1. 4. -1. 0 -2. 0 1. | MAT = | 3. 3. -1. 4. -1. 0 0 1. | | 0 0 0 -1. 4. -1. 3. 1. | | 0 0 -2. 0 -1. 4. 0 1. | | 5. 0 0 0 3. 0 4. 1. | | 0 0 1. 1. 1. 1. 1. 4. |MAT = A + D + A(T), where D is the main diagonal and A(T) is the transposed matrix of A. To store this symmetrical matrix, we can use the storage scheme of example 4.1. All we have to do is to set lsym =.true.
Program by L. Grosz, H. Haefner, C. Roll, P. Sternecker, R. Weiss 1989-96. Please contact L.Grosz.
By L. Grosz, Auckland, 4 June 2000.