/[escript]/branches/doubleplusgood/paso/src/SchurComplement.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/SchurComplement.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26