source: issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp@ 23641

Last change on this file since 23641 was 23641, checked in by Mathieu Morlighem, 6 years ago

CHG: speeding up allocation without MPI_Allreduce
File size: 9.0 KB
Line 
1/*!\file AllocateSystemMatricesx
2 * \brief retrieve vector from inputs in elements
3 */
4
5#include "./AllocateSystemMatricesx.h"
6#include "../../shared/shared.h"
7#include "../../toolkits/toolkits.h"
8
9void AllocateSystemMatricesx(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf,FemModel* femmodel){
10
11 /*Intermediary*/
12 int fsize,ssize,flocalsize,slocalsize;
13 int connectivity, numberofdofspernode;
14 int m,n,M,N;
15 int *d_nnz = NULL;
16 int *o_nnz = NULL;
17
18 /*output*/
19 Matrix<IssmDouble> *Kff = NULL;
20 Matrix<IssmDouble> *Kfs = NULL;
21 Vector<IssmDouble> *pf = NULL;
22 Vector<IssmDouble> *df = NULL;
23
24 bool oldalloc=false;
25 char* toolkittype=NULL;
26
27 /*retrieve parameters: */
28 femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
29
30 /*retrieve node info*/
31 fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
32 ssize = femmodel->nodes->NumberOfDofs(SsetEnum);
33 flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
34 slocalsize = femmodel->nodes->NumberOfDofsLocal(SsetEnum);
35
36 numberofdofspernode=femmodel->nodes->MaxNumDofs(GsetEnum);
37
38 /*if our matrices are coming from issm, we don't do dynamic allocation like Petsc
39 * does, and this routine is essentially useless. Force standard alloc in this case: */
40 toolkittype=ToolkitOptions::GetToolkitType();
41
42 if(oldalloc){
43 if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
44 if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
45 if(pdf) df =new Vector<IssmDouble>(fsize);
46 if(ppf) pf =new Vector<IssmDouble>(fsize);
47 }
48 else{
49 if(pKff){
50 m=flocalsize; n=flocalsize; /*local sizes*/
51 M=fsize; N=fsize; /*global sizes*/
52 if(strcmp(toolkittype,"issm")==0){
53 Kff=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
54 }
55 else{
56 MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,FsetEnum);
57 Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
58 xDelete<int>(d_nnz);
59 xDelete<int>(o_nnz);
60 }
61 }
62 if(pKfs){
63 m=flocalsize; n=slocalsize; /*local sizes*/
64 M=fsize; N=ssize; /*global sizes*/
65 if(strcmp(toolkittype,"issm")==0){
66 Kfs=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
67 }
68 else{
69 MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,SsetEnum);
70 Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
71 xDelete<int>(d_nnz);
72 xDelete<int>(o_nnz);
73 }
74 }
75 if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
76 if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
77 }
78
79 /*Free ressources: */
80 xDelete<char>(toolkittype);
81
82 /*Allocate output pointers*/
83 if(pKff) *pKff = Kff;
84 if(pKfs) *pKfs = Kfs;
85 if(pdf) *pdf = df;
86 if(ppf) *ppf = pf;
87}
88
89void MatrixNonzeros(int** pd_nnz,int** po_nnz,FemModel* femmodel,int set1enum,int set2enum){
90
91 /*Intermediary*/
92 int i,j,k,index,offset,count;
93 int d_nz,o_nz;
94 Element *element = NULL;
95 Load *load = NULL;
96 int *head_e = NULL;
97 int *next_e = NULL;
98 int *count2offset_e = NULL;
99 int *head_l = NULL;
100 int *next_l = NULL;
101 int *count2offset_l = NULL;
102 int *lidlist = NULL;
103
104 /*output*/
105 int *d_nnz = NULL;
106 int *o_nnz = NULL;
107
108 /*Get vector size and number of nodes*/
109 int numnodes = femmodel->nodes->NumberOfNodes();
110 int localmasters = femmodel->nodes->NumberOfNodesLocal();
111 int localnumnodes = femmodel->nodes->Size();
112 int numberofdofspernode = femmodel->nodes->MaxNumDofs(GsetEnum);
113 int M = femmodel->nodes->NumberOfDofs(set1enum);
114 int N = femmodel->nodes->NumberOfDofs(set2enum);
115 int m = femmodel->nodes->NumberOfDofsLocal(set1enum);
116 int n = femmodel->nodes->NumberOfDofsLocal(set2enum);
117 int numnodesperelement = femmodel->elements->MaxNumNodes();
118 int numnodesperload = femmodel->loads->MaxNumNodes();
119
120 /*First, we are building chaining vectors so that we know what nodes are
121 * connected to what elements. These vectors are such that:
122 * for(int i=head[id];i!=-1;i=next[i])
123 * will loop over all the elements that are connected to the node number
124 * id*/
125 head_e = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
126 next_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
127 count2offset_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
128
129 k=0;
130 for(i=0;i<femmodel->elements->Size();i++){
131 element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
132 lidlist = xNew<int>(element->GetNumberOfNodes());
133 element->GetNodesLidList(lidlist);
134
135 for(j=0;j<element->GetNumberOfNodes();j++){
136 index = lidlist[j];
137 _assert_(index>=0 && index<numnodes);
138
139 count2offset_e[k]=i;
140 next_e[k]=head_e[index];
141 head_e[index]=k++;
142 }
143 for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
144
145 xDelete<int>(lidlist);
146 }
147
148 /*Chain for loads*/
149 head_l = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
150 next_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
151 count2offset_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
152 k=0;
153 for(i=0;i<femmodel->loads->Size();i++){
154 load = xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
155 lidlist = xNew<int>(load->GetNumberOfNodes());
156 load->GetNodesLidList(lidlist);
157
158 for(j=0;j<load->GetNumberOfNodes();j++){
159 index = lidlist[j];
160 _assert_(index>=0 && index<numnodes);
161
162 count2offset_l[k]=i;
163 next_l[k]=head_l[index];
164 head_l[index]=k++;
165 }
166 for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
167
168 xDelete<int>(lidlist);
169 }
170
171 /*OK now count number of dofs and flag each nodes for each node i*/
172 bool *flags = xNew<bool>(localnumnodes);
173 int *flagsindices = xNew<int>(localnumnodes);
174 int *d_connectivity = xNewZeroInit<int>(localnumnodes);
175 int *o_connectivity = xNewZeroInit<int>(localnumnodes);
176 int *all_connectivity_clone = xNewZeroInit<int>(numnodes);
177
178 Vector<IssmDouble> *connectivity_clone= new Vector<IssmDouble>(localmasters,numnodes);
179
180 /*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
181 * to reset in flagsindices, initialized with -1*/
182 for(i = 0;i<localnumnodes;i++) flags[i] = false;
183 for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
184
185 /*Create connectivity vector*/
186 for(i=0;i<femmodel->nodes->Size();i++){
187 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
188
189 /*Reinitialize flags to false*/
190 j=0;
191 while(j<localnumnodes){
192 if(flagsindices[j]>=0){
193 flags[flagsindices[j]] = false;
194 flagsindices[j] = -1;
195 j++;
196 }
197 else{
198 break;
199 }
200 }
201
202 //for(j=0;j<localnumnodes;j++) flags[j]=false;
203
204 /*Loop over elements that hold node number i*/
205 //if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
206 // printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
207 //}
208 for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
209 offset=count2offset_e[j];
210 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(offset));
211 element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
212 if(node->IsClone()){
213 connectivity_clone->SetValue(node->Pid(),d_nz+o_nz,ADD_VAL);
214 }
215 else{
216 d_connectivity[node->Lid()]+=d_nz;
217 o_connectivity[node->Lid()]+=o_nz;
218 }
219 }
220 for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
221 offset=count2offset_l[j];
222 load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(offset));
223 load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
224 if(node->IsClone()){
225 connectivity_clone->SetValue(node->Pid(),d_nz+o_nz,ADD_VAL);
226 }
227 else{
228 d_connectivity[node->Lid()]+=d_nz;
229 o_connectivity[node->Lid()]+=o_nz;
230 }
231 }
232 }
233 xDelete<bool>(flags);
234 xDelete<int>(flagsindices);
235 xDelete<int>(count2offset_e);
236 xDelete<int>(head_e);
237 xDelete<int>(next_e);
238 xDelete<int>(count2offset_l);
239 xDelete<int>(head_l);
240 xDelete<int>(next_l);
241
242 /*sum over all cpus*/
243 connectivity_clone->Assemble();
244 IssmDouble* serial_connectivity_clone=NULL;
245 femmodel->GetLocalVectorWithClonesVertices(&serial_connectivity_clone,connectivity_clone);
246 delete connectivity_clone;
247
248 if(set1enum==FsetEnum){
249 count=0;
250 d_nnz=xNew<int>(m);
251 o_nnz=xNew<int>(m);
252 for(i=0;i<femmodel->nodes->Size();i++){
253 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
254 if(!node->IsClone()){
255 for(j=0;j<node->fsize;j++){
256 _assert_(count<m);
257 d_nnz[count]=numberofdofspernode*(d_connectivity[node->Lid()] + reCast<int>(serial_connectivity_clone[node->Lid()]));
258 o_nnz[count]=numberofdofspernode*(o_connectivity[node->Lid()] + reCast<int>(serial_connectivity_clone[node->Lid()]));
259 if(d_nnz[count]>n) d_nnz[count]=n;
260 if(o_nnz[count]>N-n) o_nnz[count]=N-n;
261 count++;
262 }
263 }
264 }
265 _assert_(m==count);
266 }
267 else{
268 _error_("STOP not implemented");
269 }
270 xDelete<int>(d_connectivity);
271 xDelete<int>(o_connectivity);
272 xDelete<IssmDouble>(serial_connectivity_clone);
273
274 /*Allocate ouptput pointer*/
275 *pd_nnz=d_nnz;
276 *po_nnz=o_nnz;
277}
Note: See TracBrowser for help on using the repository browser.