Annotation of /trunk/paso/src/Solver_SchurComplement.c

Revision 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (11 years, 9 months ago) by phornby
Original Path: temp_trunk_copy/paso/src/Solver_SchurComplement.c
File MIME type: text/plain
File size: 10657 byte(s)
```Make a temp copy of the trunk before checking in the windows changes

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

Properties

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