# Contents of /trunk/paso/src/Solver_SchurComplement.c

Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/paso/src/Solvers/Solver_SchurComplement.c
File MIME type: text/plain
File size: 10213 byte(s)
```move all directories from trunk/esys2 into trunk and remove esys2

```
 1 /* \$Id\$ */ 2 3 /**************************************************************/ 4 5 /* Paso: updates A_CC <- ACC-ACF AFF^{-1} AFC */ 6 7 /* no check of consistency of matrices !!!! */ 8 9 /**************************************************************/ 10 11 /* Copyrights by ACcESS Australia 2003,2004,2005 */ 12 /* Author: gross@access.edu.au */ 13 14 /**************************************************************/ 15 16 #include "Paso.h" 17 #include "SystemMatrix.h" 18 #include "Solver.h" 19 20 /**************************************************************/ 21 22 23 24 void Paso_Solver_updateIncompleteSchurComplement(Paso_SystemMatrix* A_CC,Paso_SystemMatrix *A_CF,double* invA_FF,index_t* A_FF_pivot,Paso_SystemMatrix *A_FC) { 25 index_t iPtr_CC,*index_CC,col_CF,col_FC, *where_p,iPtr_CC_2,i,iPtr_CF,iPtr_FC; 26 dim_t index_CC_len; 27 bool_t set_A; 28 dim_t n_rows=A_CC->num_rows; 29 dim_t n_block=A_CC->row_block_size; 30 register double A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33, 31 invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33, 32 A11,A21,A31,A12,A22,A32,A13,A23,A33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33; 33 if (n_block==1) { 34 #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,col_CF,set_A,iPtr_CF,iPtr_FC,col_FC,where_p,A11) schedule(static) 35 for (i = 0; i < n_rows;++i) { 36 iPtr_CC=A_CC->pattern->ptr[i]; 37 index_CC=&(A_CC->pattern->index[iPtr_CC]); 38 index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]); 39 /* now we run through the columns of A_CF in row i */ 40 for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) { 41 col_CF=A_CF->pattern->index[iPtr_CF]; 42 set_A=TRUE; 43 for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) { 44 col_FC=A_FC->pattern->index[iPtr_FC]; 45 /* is (i,col_FC) in the shape of A_CC ? */ 46 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex); 47 if (where_p!=NULL) { 48 if (set_A) { 49 A11=A_CF->val[iPtr_CF]*invA_FF[col_CF]; 50 set_A=FALSE; 51 } 52 A_CC->val[iPtr_CC+(index_t)(where_p-index_CC)]-=A11*A_FC->val[iPtr_FC]; 53 } 54 } /* end of iPtr_FC loop */ 55 } /* end of iPtr_CF loop */ 56 } /* end of irow loop */ 57 } else if (n_block==2) { 58 #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_12,A_CF_22,invA_FF_11,invA_FF_21,invA_FF_12,invA_FF_22,A11,A21,A12,A22,A_FC_11,A_FC_21,A_FC_12,A_FC_22) schedule(static) 59 for (i = 0; i < n_rows;++i) { 60 iPtr_CC=A_CC->pattern->ptr[i]; 61 index_CC=&(A_CC->pattern->index[iPtr_CC]); 62 index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]); 63 /* now we run through the columns of A_CF in row i */ 64 for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) { 65 col_CF=A_CF->pattern->index[iPtr_CF]; 66 set_A=TRUE; 67 for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) { 68 col_FC=A_FC->pattern->index[iPtr_FC]; 69 /* is (i,col_FC) in the shape of A_CC ? */ 70 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex); 71 if (where_p!=NULL) { 72 iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC); 73 /* this calculutes A_CF*invA_FF(i,col_CF) */ 74 if (set_A) { 75 A_CF_11=A_CF->val[4*iPtr_CF ]; 76 A_CF_21=A_CF->val[4*iPtr_CF+1]; 77 A_CF_12=A_CF->val[4*iPtr_CF+2]; 78 A_CF_22=A_CF->val[4*iPtr_CF+3]; 79 80 invA_FF_11=invA_FF[4*col_CF ]; 81 invA_FF_21=invA_FF[4*col_CF+1]; 82 invA_FF_12=invA_FF[4*col_CF+2]; 83 invA_FF_22=invA_FF[4*col_CF+3]; 84 85 A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21; 86 A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21; 87 A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22; 88 A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22; 89 90 set_A=FALSE; 91 } 92 93 A_FC_11=A_FC->val[4*iPtr_FC ]; 94 A_FC_21=A_FC->val[4*iPtr_FC+1]; 95 A_FC_12=A_FC->val[4*iPtr_FC+2]; 96 A_FC_22=A_FC->val[4*iPtr_FC+3]; 97 98 A_CC->val[4*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21; 99 A_CC->val[4*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21; 100 A_CC->val[4*iPtr_CC_2+2]-=A11*A_FC_12+A12*A_FC_22; 101 A_CC->val[4*iPtr_CC_2+3]-=A21*A_FC_12+A22*A_FC_22; 102 103 } 104 } /* end of iPtr_FC loop */ 105 } /* end of iPtr_CF loop */ 106 } /* end of irow loop */ 107 } else if (n_block==3) { 108 #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,A11,A21,A31,A12,A22,A32,A13,A23,A33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33) schedule(static) 109 for (i = 0; i < n_rows;++i) { 110 iPtr_CC=A_CC->pattern->ptr[i]; 111 index_CC=&(A_CC->pattern->index[iPtr_CC]); 112 index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]); 113 /* now we run through the columns of A_CF in row i */ 114 for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) { 115 col_CF=A_CF->pattern->index[iPtr_CF]; 116 set_A=TRUE; 117 for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) { 118 col_FC=A_FC->pattern->index[iPtr_FC]; 119 /* is (i,col_FC) in the shape of A_CC ? */ 120 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex); 121 if (where_p!=NULL) { 122 iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC); 123 /* this calculutes A_CF*invA_FF(i,col_CF) */ 124 if (set_A) { 125 A_CF_11=A_CF->val[9*iPtr_CF ]; 126 A_CF_21=A_CF->val[9*iPtr_CF+1]; 127 A_CF_31=A_CF->val[9*iPtr_CF+2]; 128 A_CF_12=A_CF->val[9*iPtr_CF+3]; 129 A_CF_22=A_CF->val[9*iPtr_CF+4]; 130 A_CF_32=A_CF->val[9*iPtr_CF+5]; 131 A_CF_13=A_CF->val[9*iPtr_CF+6]; 132 A_CF_23=A_CF->val[9*iPtr_CF+7]; 133 A_CF_33=A_CF->val[9*iPtr_CF+8]; 134 135 invA_FF_11=invA_FF[9*col_CF ]; 136 invA_FF_21=invA_FF[9*col_CF+1]; 137 invA_FF_31=invA_FF[9*col_CF+2]; 138 invA_FF_12=invA_FF[9*col_CF+3]; 139 invA_FF_22=invA_FF[9*col_CF+4]; 140 invA_FF_32=invA_FF[9*col_CF+5]; 141 invA_FF_13=invA_FF[9*col_CF+6]; 142 invA_FF_23=invA_FF[9*col_CF+7]; 143 invA_FF_33=invA_FF[9*col_CF+8]; 144 145 A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21+A_CF_13*invA_FF_31; 146 A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21+A_CF_23*invA_FF_31; 147 A31=A_CF_31*invA_FF_11+A_CF_32*invA_FF_21+A_CF_33*invA_FF_31; 148 A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22+A_CF_13*invA_FF_32; 149 A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22+A_CF_23*invA_FF_32; 150 A32=A_CF_31*invA_FF_12+A_CF_32*invA_FF_22+A_CF_33*invA_FF_32; 151 A13=A_CF_11*invA_FF_13+A_CF_12*invA_FF_23+A_CF_13*invA_FF_33; 152 A23=A_CF_21*invA_FF_13+A_CF_22*invA_FF_23+A_CF_23*invA_FF_33; 153 A33=A_CF_31*invA_FF_13+A_CF_32*invA_FF_23+A_CF_33*invA_FF_33; 154 155 set_A=FALSE; 156 } 157 158 A_FC_11=A_FC->val[9*iPtr_FC ]; 159 A_FC_21=A_FC->val[9*iPtr_FC+1]; 160 A_FC_31=A_FC->val[9*iPtr_FC+2]; 161 A_FC_12=A_FC->val[9*iPtr_FC+3]; 162 A_FC_22=A_FC->val[9*iPtr_FC+4]; 163 A_FC_32=A_FC->val[9*iPtr_FC+5]; 164 A_FC_13=A_FC->val[9*iPtr_FC+6]; 165 A_FC_23=A_FC->val[9*iPtr_FC+7]; 166 A_FC_33=A_FC->val[9*iPtr_FC+8]; 167 168 A_CC->val[9*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21+A13*A_FC_31; 169 A_CC->val[9*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21+A23*A_FC_31; 170 A_CC->val[9*iPtr_CC_2+2]-=A31*A_FC_11+A32*A_FC_21+A33*A_FC_31; 171 A_CC->val[9*iPtr_CC_2+3]-=A11*A_FC_12+A12*A_FC_22+A13*A_FC_32; 172 A_CC->val[9*iPtr_CC_2+4]-=A21*A_FC_12+A22*A_FC_22+A23*A_FC_32; 173 A_CC->val[9*iPtr_CC_2+5]-=A31*A_FC_12+A32*A_FC_22+A33*A_FC_32; 174 A_CC->val[9*iPtr_CC_2+6]-=A11*A_FC_13+A12*A_FC_23+A13*A_FC_33; 175 A_CC->val[9*iPtr_CC_2+7]-=A21*A_FC_13+A22*A_FC_23+A23*A_FC_33; 176 A_CC->val[9*iPtr_CC_2+8]-=A31*A_FC_13+A32*A_FC_23+A33*A_FC_33; 177 178 } 179 } /* end of iPtr_FC loop */ 180 } /* end of iPtr_CF loop */ 181 } /* end of irow loop */ 182 } 183 } 184 /* 185 * \$Log\$ 186 * Revision 1.2 2005/09/15 03:44:40 jgs 187 * Merge of development branch dev-02 back to main trunk on 2005-09-15 188 * 189 * Revision 1.1.2.1 2005/09/05 06:29:50 gross 190 * These files have been extracted from finley to define a stand alone libray for iterative 191 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but 192 * has not been tested yet. 193 * 194 * 195 */

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision