/[escript]/trunk/paso/src/SchurComplement.cpp
ViewVC logotype

Annotation of /trunk/paso/src/SchurComplement.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (hide annotations)
Mon Jan 20 03:37:18 2020 UTC (2 months, 1 week ago) by uqaeller
File size: 10195 byte(s)
Updated the copyright header.


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/paso/src/SchurComplement.cpp:5567-5588 /branches/amg_from_3530/paso/src/SchurComplement.cpp:3531-3826 /branches/complex/paso/src/SchurComplement.cpp:5866-5937 /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 /branches/trilinos_from_5897/paso/src/SchurComplement.cpp:5898-6118 /release/3.0/paso/src/SchurComplement.cpp:2591-2601 /release/4.0/paso/src/SchurComplement.cpp:5380-5406 /trunk/paso/src/SchurComplement.cpp:4257-4344 /trunk/ripley/test/python/paso/src/SchurComplement.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26