# Contents of /trunk/paso/src/SchurComplement.cpp

Revision 4873 - (show annotations)
Wed Apr 16 06:38:51 2014 UTC (5 years, 5 months ago) by caltinay
File size: 10559 byte(s)
whitespace only changes.

 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2014 by University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 * 15 *****************************************************************************/ 16 17 18 /****************************************************************************/ 19 20 /* Paso: updates A_CC <- ACC-ACF AFF^{-1} AFC */ 21 22 /* no check of consistency of matrices !!!! */ 23 24 /****************************************************************************/ 25 26 /* Copyrights by ACcESS Australia 2003,2004,2005 */ 27 /* Author: Lutz Gross, l.gross@uq.edu.au */ 28 29 /****************************************************************************/ 30 31 #include "Paso.h" 32 #include "SparseMatrix.h" 33 #include "Solver.h" 34 35 /****************************************************************************/ 36 37 namespace paso { 38 39 void Solver_updateIncompleteSchurComplement(SparseMatrix_ptr A_CC, 40 SparseMatrix_ptr A_CF, double* invA_FF, index_t* A_FF_pivot, 41 SparseMatrix_ptr A_FC) 42 { 43 index_t iPtr_CC,*index_CC,col_CF,col_FC, *where_p,iPtr_CC_2,i,iPtr_CF,iPtr_FC; 44 dim_t index_CC_len; 45 bool set_A; 46 dim_t n_loc_rows=A_CC->numRows; 47 dim_t n_block=A_CC->row_block_size; 48 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, 49 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, 50 A11=0,A21=0,A31=0,A12=0,A22=0,A32=0,A13=0,A23=0,A33=0,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; 51 if (n_block==1) { 52 #pragma omp parallel for firstprivate(A11) private(i,iPtr_CC,index_CC,index_CC_len,col_CF,set_A,iPtr_CF,iPtr_FC,col_FC,where_p) schedule(static) 53 for (i = 0; i < n_loc_rows;++i) { 54 iPtr_CC=A_CC->pattern->ptr[i]; 55 index_CC=&(A_CC->pattern->index[iPtr_CC]); 56 index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]); 57 /* now we run through the columns of A_CF in row i */ 58 for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) { 59 col_CF=A_CF->pattern->index[iPtr_CF]; 60 set_A=true; 61 for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) { 62 col_FC=A_FC->pattern->index[iPtr_FC]; 63 /* is (i,col_FC) in the shape of A_CC ? */ 64 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),util::comparIndex); 65 if (where_p!=NULL) { 66 if (set_A) { 67 A11=A_CF->val[iPtr_CF]*invA_FF[col_CF]; 68 set_A=false; 69 } 70 A_CC->val[iPtr_CC+(index_t)(where_p-index_CC)]-=A11*A_FC->val[iPtr_FC]; 71 } 72 } /* end of iPtr_FC loop */ 73 } /* end of iPtr_CF loop */ 74 } /* end of irow loop */ 75 } else if (n_block==2) { 76 #pragma omp parallel for firstprivate(A11,A21,A12,A22) 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,A_FC_11,A_FC_21,A_FC_12,A_FC_22) schedule(static) 77 for (i = 0; i < n_loc_rows;++i) { 78 iPtr_CC=A_CC->pattern->ptr[i]; 79 index_CC=&(A_CC->pattern->index[iPtr_CC]); 80 index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]); 81 /* now we run through the columns of A_CF in row i */ 82 for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) { 83 col_CF=A_CF->pattern->index[iPtr_CF]; 84 set_A=true; 85 for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) { 86 col_FC=A_FC->pattern->index[iPtr_FC]; 87 /* is (i,col_FC) in the shape of A_CC ? */ 88 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),util::comparIndex); 89 if (where_p!=NULL) { 90 iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC); 91 /* this calculates A_CF*invA_FF(i,col_CF) */ 92 if (set_A) { 93 A_CF_11=A_CF->val[4*iPtr_CF ]; 94 A_CF_21=A_CF->val[4*iPtr_CF+1]; 95 A_CF_12=A_CF->val[4*iPtr_CF+2]; 96 A_CF_22=A_CF->val[4*iPtr_CF+3]; 97 98 invA_FF_11=invA_FF[4*col_CF ]; 99 invA_FF_21=invA_FF[4*col_CF+1]; 100 invA_FF_12=invA_FF[4*col_CF+2]; 101 invA_FF_22=invA_FF[4*col_CF+3]; 102 103 A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21; 104 A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21; 105 A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22; 106 A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22; 107 108 set_A=false; 109 } 110 111 A_FC_11=A_FC->val[4*iPtr_FC ]; 112 A_FC_21=A_FC->val[4*iPtr_FC+1]; 113 A_FC_12=A_FC->val[4*iPtr_FC+2]; 114 A_FC_22=A_FC->val[4*iPtr_FC+3]; 115 116 A_CC->val[4*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21; 117 A_CC->val[4*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21; 118 A_CC->val[4*iPtr_CC_2+2]-=A11*A_FC_12+A12*A_FC_22; 119 A_CC->val[4*iPtr_CC_2+3]-=A21*A_FC_12+A22*A_FC_22; 120 121 } 122 } /* end of iPtr_FC loop */ 123 } /* end of iPtr_CF loop */ 124 } /* end of irow loop */ 125 } else if (n_block==3) { 126 #pragma omp parallel for firstprivate(A11,A21,A31,A12,A22,A32,A13,A23,A33) 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,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) 127 for (i = 0; i < n_loc_rows;++i) { 128 iPtr_CC=A_CC->pattern->ptr[i]; 129 index_CC=&(A_CC->pattern->index[iPtr_CC]); 130 index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]); 131 /* now we run through the columns of A_CF in row i */ 132 for (iPtr_CF = A_CF->pattern->ptr[i]; iPtr_CF < A_CF->pattern->ptr[i + 1]; ++iPtr_CF) { 133 col_CF=A_CF->pattern->index[iPtr_CF]; 134 set_A = true; 135 for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) { 136 col_FC=A_FC->pattern->index[iPtr_FC]; 137 /* is (i,col_FC) in the shape of A_CC ? */ 138 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),util::comparIndex); 139 if (where_p!=NULL) { 140 iPtr_CC_2=iPtr_CC+(index_t)(where_p-index_CC); 141 /* this calculates A_CF*invA_FF(i,col_CF) */ 142 if (set_A) { 143 A_CF_11=A_CF->val[9*iPtr_CF ]; 144 A_CF_21=A_CF->val[9*iPtr_CF+1]; 145 A_CF_31=A_CF->val[9*iPtr_CF+2]; 146 A_CF_12=A_CF->val[9*iPtr_CF+3]; 147 A_CF_22=A_CF->val[9*iPtr_CF+4]; 148 A_CF_32=A_CF->val[9*iPtr_CF+5]; 149 A_CF_13=A_CF->val[9*iPtr_CF+6]; 150 A_CF_23=A_CF->val[9*iPtr_CF+7]; 151 A_CF_33=A_CF->val[9*iPtr_CF+8]; 152 153 invA_FF_11=invA_FF[9*col_CF ]; 154 invA_FF_21=invA_FF[9*col_CF+1]; 155 invA_FF_31=invA_FF[9*col_CF+2]; 156 invA_FF_12=invA_FF[9*col_CF+3]; 157 invA_FF_22=invA_FF[9*col_CF+4]; 158 invA_FF_32=invA_FF[9*col_CF+5]; 159 invA_FF_13=invA_FF[9*col_CF+6]; 160 invA_FF_23=invA_FF[9*col_CF+7]; 161 invA_FF_33=invA_FF[9*col_CF+8]; 162 163 A11=A_CF_11*invA_FF_11+A_CF_12*invA_FF_21+A_CF_13*invA_FF_31; 164 A21=A_CF_21*invA_FF_11+A_CF_22*invA_FF_21+A_CF_23*invA_FF_31; 165 A31=A_CF_31*invA_FF_11+A_CF_32*invA_FF_21+A_CF_33*invA_FF_31; 166 A12=A_CF_11*invA_FF_12+A_CF_12*invA_FF_22+A_CF_13*invA_FF_32; 167 A22=A_CF_21*invA_FF_12+A_CF_22*invA_FF_22+A_CF_23*invA_FF_32; 168 A32=A_CF_31*invA_FF_12+A_CF_32*invA_FF_22+A_CF_33*invA_FF_32; 169 A13=A_CF_11*invA_FF_13+A_CF_12*invA_FF_23+A_CF_13*invA_FF_33; 170 A23=A_CF_21*invA_FF_13+A_CF_22*invA_FF_23+A_CF_23*invA_FF_33; 171 A33=A_CF_31*invA_FF_13+A_CF_32*invA_FF_23+A_CF_33*invA_FF_33; 172 173 set_A=false; 174 } 175 176 A_FC_11=A_FC->val[9*iPtr_FC ]; 177 A_FC_21=A_FC->val[9*iPtr_FC+1]; 178 A_FC_31=A_FC->val[9*iPtr_FC+2]; 179 A_FC_12=A_FC->val[9*iPtr_FC+3]; 180 A_FC_22=A_FC->val[9*iPtr_FC+4]; 181 A_FC_32=A_FC->val[9*iPtr_FC+5]; 182 A_FC_13=A_FC->val[9*iPtr_FC+6]; 183 A_FC_23=A_FC->val[9*iPtr_FC+7]; 184 A_FC_33=A_FC->val[9*iPtr_FC+8]; 185 186 A_CC->val[9*iPtr_CC_2 ]-=A11*A_FC_11+A12*A_FC_21+A13*A_FC_31; 187 A_CC->val[9*iPtr_CC_2+1]-=A21*A_FC_11+A22*A_FC_21+A23*A_FC_31; 188 A_CC->val[9*iPtr_CC_2+2]-=A31*A_FC_11+A32*A_FC_21+A33*A_FC_31; 189 A_CC->val[9*iPtr_CC_2+3]-=A11*A_FC_12+A12*A_FC_22+A13*A_FC_32; 190 A_CC->val[9*iPtr_CC_2+4]-=A21*A_FC_12+A22*A_FC_22+A23*A_FC_32; 191 A_CC->val[9*iPtr_CC_2+5]-=A31*A_FC_12+A32*A_FC_22+A33*A_FC_32; 192 A_CC->val[9*iPtr_CC_2+6]-=A11*A_FC_13+A12*A_FC_23+A13*A_FC_33; 193 A_CC->val[9*iPtr_CC_2+7]-=A21*A_FC_13+A22*A_FC_23+A23*A_FC_33; 194 A_CC->val[9*iPtr_CC_2+8]-=A31*A_FC_13+A32*A_FC_23+A33*A_FC_33; 195 196 } 197 } /* end of iPtr_FC loop */ 198 } /* end of iPtr_CF loop */ 199 } /* end of irow loop */ 200 } 201 } 202 203 } // namespace paso 204

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/amg_from_3530/paso/src/SchurComplement.cpp:3531-3826 /branches/lapack2681/paso/src/SchurComplement.cpp:2682-2741 /branches/pasowrap/paso/src/SchurComplement.cpp:3661-3674 /branches/py3_attempt2/paso/src/SchurComplement.cpp:3871-3891 /branches/restext/paso/src/SchurComplement.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/SchurComplement.cpp:3669-3791 /branches/stage3.0/paso/src/SchurComplement.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/SchurComplement.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/SchurComplement.cpp:3517-3974 /release/3.0/paso/src/SchurComplement.cpp:2591-2601 /trunk/paso/src/SchurComplement.cpp:4257-4344 /trunk/ripley/test/python/paso/src/SchurComplement.cpp:3480-3515