/[escript]/trunk-mpi-branch/finley/src/Mesh_getPattern.c
ViewVC logotype

Contents of /trunk-mpi-branch/finley/src/Mesh_getPattern.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1096 - (show annotations)
Mon Apr 16 22:59:33 2007 UTC (14 years ago) by ksteube
File MIME type: text/plain
File size: 8008 byte(s)
MPI implicit solver example run_simplesolve.py now compiling and
running successfully on one CPU of ess.

Adjusted SConscript, removed some debug print statements and removed
some partially implemented TRILINOS calls.


1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13 /**************************************************************/
14
15 /* Finley: Mesh */
16
17 /**************************************************************/
18
19 /* Author: gross@access.edu.au */
20 /* Version: $Id$ */
21
22 /**************************************************************/
23
24 #include "Mesh.h"
25 #include "IndexList.h"
26
27 /**************************************************************/
28
29 /* returns a reference to the matrix pattern */
30
31 Paso_SystemMatrixPattern* Finley_getPattern(Finley_Mesh *mesh,bool_t reduce_row_order, bool_t reduce_col_order) {
32 Paso_SystemMatrixPattern *out=NULL;
33 Finley_resetError();
34 /* make sure that the requested pattern is available */
35 if (reduce_row_order) {
36 if (reduce_col_order) {
37 if (mesh->ReducedReducedPattern==NULL) mesh->ReducedReducedPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
38 } else {
39 if (mesh->ReducedFullPattern==NULL) mesh->ReducedFullPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
40 }
41 } else {
42 if (reduce_col_order) {
43 if (mesh->FullReducedPattern==NULL) mesh->FullReducedPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
44 } else {
45 if (mesh->FullFullPattern==NULL) mesh->FullFullPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
46 }
47 }
48 if (Finley_noError()) {
49 if (reduce_row_order) {
50 if (reduce_col_order) {
51 out=Paso_SystemMatrixPattern_reference(mesh->ReducedReducedPattern);
52 } else {
53 out=Paso_SystemMatrixPattern_reference(mesh->ReducedFullPattern);
54 }
55 } else {
56 if (reduce_col_order) {
57 out=Paso_SystemMatrixPattern_reference(mesh->FullReducedPattern);
58 } else {
59 out=Paso_SystemMatrixPattern_reference(mesh->FullFullPattern);
60 }
61 }
62 }
63 return out;
64 }
65 Paso_SystemMatrixPattern* Finley_makePattern(Finley_Mesh *mesh,bool_t reduce_row_order, bool_t reduce_col_order) {
66 double time0;
67 Paso_SystemMatrixPattern *pattern = NULL;
68 dim_t i,n, numHops;
69 index_t s,itmp,*rowLabel=NULL,*colLabel=NULL, *ptr=NULL, *index=NULL, *hop=NULL;
70 Finley_IndexList* index_list=NULL;
71 Finley_resetError();
72 Finley_NodeDistribution *row_degreeOfFreedomDistribution=NULL;
73 Finley_NodeDistribution *col_degreeOfFreedomDistribution=NULL;
74 Paso_Distribution *row_distribution=NULL;
75 Paso_Distribution *col_distribution=NULL;
76
77 time0=Finley_timer();
78
79 #ifdef PASO_MPI
80 if (reduce_col_order) {
81 n=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
82 colLabel=mesh->Nodes->reducedDegreeOfFreedom;
83 } else {
84 n=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
85 colLabel=mesh->Nodes->degreeOfFreedom;
86 }
87
88 if (reduce_row_order) {
89 n=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
90 rowLabel=mesh->Nodes->reducedDegreeOfFreedom;
91 } else {
92 n=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
93 rowLabel=mesh->Nodes->degreeOfFreedom;
94 }
95 #else
96 if (reduce_col_order) {
97 n=mesh->Nodes->reducedNumDegreesOfFreedom;
98 colLabel=mesh->Nodes->reducedDegreeOfFreedom;
99 } else {
100 n=mesh->Nodes->numDegreesOfFreedom;
101 colLabel=mesh->Nodes->degreeOfFreedom;
102 }
103
104 if (reduce_row_order) {
105 n=mesh->Nodes->reducedNumDegreesOfFreedom;
106 rowLabel=mesh->Nodes->reducedDegreeOfFreedom;
107 } else {
108 n=mesh->Nodes->numDegreesOfFreedom;
109 rowLabel=mesh->Nodes->degreeOfFreedom;
110 }
111 #endif
112
113 index_list=TMPMEMALLOC(n,Finley_IndexList);
114 ptr=MEMALLOC(n+1,index_t);
115 if (! (Finley_checkPtr(index_list) || Finley_checkPtr(ptr)) ) {
116
117 #pragma omp parallel private(i,s,itmp)
118 {
119 #pragma omp for schedule(static)
120 for(i=0;i<n;++i) {
121 index_list[i].extension=NULL;
122 index_list[i].n=0;
123 }
124 /* insert contributions from element matrices into colums index index_list: */
125 Finley_IndexList_insertElements(index_list,mesh,mesh->Elements,
126 reduce_row_order,rowLabel,reduce_col_order,colLabel);
127 Finley_IndexList_insertElements(index_list,mesh,mesh->FaceElements,
128 reduce_row_order,rowLabel,reduce_col_order,colLabel);
129 Finley_IndexList_insertElements(index_list,mesh,mesh->ContactElements,
130 reduce_row_order,rowLabel,reduce_col_order,colLabel);
131 Finley_IndexList_insertElements(index_list,mesh,mesh->Points,
132 reduce_row_order,rowLabel,reduce_col_order,colLabel);
133 /* get the number of connections per row */
134 #pragma omp for schedule(static)
135 for(i=0;i<n;++i) {
136 ptr[i]=Finley_IndexList_count(&index_list[i]);
137 }
138 /* accumalate ptr */
139 /* OMP */
140 #pragma omp master
141 {
142 s=0;
143 for(i=0;i<n;++i) {
144 itmp=ptr[i];
145 ptr[i]=s;
146 s+=itmp;
147 }
148 ptr[n]=s;
149 /* fill index */
150 index=MEMALLOC(ptr[n],index_t);
151 }
152 #pragma omp barrier
153 if (! Finley_checkPtr(index)) {
154 #pragma omp for schedule(static)
155 for(i=0;i<n;++i) Finley_IndexList_toArray(&index_list[i],&index[ptr[i]]);
156 }
157 }
158 }
159 /* all done */
160 /* clean up */
161 if (index_list!=NULL) {
162 #pragma omp parallel for private(i)
163 for(i=0;i<n;++i) Finley_IndexList_free(index_list[i].extension);
164 }
165 TMPMEMFREE(index_list);
166 #ifdef Finley_TRACE
167 printf("timing: mesh to matrix pattern: %.4e sec\n",Finley_timer()-time0);
168 #endif
169 if (Finley_noError()) {
170 // Use a getReference method to get the DOF distributions to avoid memory leaks
171 if (reduce_col_order) {
172 col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
173 } else {
174 col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
175 }
176 if (reduce_row_order) {
177 row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
178 } else {
179 row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
180 }
181 row_distribution=Paso_Distribution_alloc(mesh->MPIInfo, row_degreeOfFreedomDistribution->vtxdist,1,0);
182 col_distribution=Paso_Distribution_alloc(mesh->MPIInfo, col_degreeOfFreedomDistribution->vtxdist,1,0);
183
184 if (Finley_noError()) {
185 Paso_SystemMatrixPattern_makeHops(PATTERN_FORMAT_DEFAULT,col_distribution,ptr,index,&numHops,&hop);
186 if (Finley_noError()) {
187 pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,
188 row_distribution,
189 col_distribution,
190 ptr,index,numHops, hop);
191 if (Finley_noError()) {
192 #ifdef PASO_MPI
193 pattern->output_node_distribution=Finley_NodeDistribution_getReference(row_degreeOfFreedomDistribution);
194 pattern->input_node_distribution=Finley_NodeDistribution_getReference(col_degreeOfFreedomDistribution);
195 #endif
196 }
197 }
198 }
199 Paso_Distribution_free(row_distribution);
200 Paso_Distribution_free(col_distribution);
201 }
202 MEMFREE(hop);
203 if (! Finley_noError()) {
204 MEMFREE(ptr);
205 MEMFREE(index);
206 }
207 return pattern;
208 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26