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

Contents of /trunk/finley/src/Mesh_getPattern.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 969 - (show annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 11 months ago) by ksteube
File MIME type: text/plain
File size: 8435 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


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 out->MPIInfo = Paso_MPIInfo_getReference( mesh->MPIInfo );
64 return out;
65 }
66 Paso_SystemMatrixPattern* Finley_makePattern(Finley_Mesh *mesh,bool_t reduce_row_order, bool_t reduce_col_order) {
67 double time0;
68 dim_t i,n;
69 index_t s,itmp,*rowLabel=NULL,*colLabel=NULL, *ptr=NULL, *index=NULL;
70 Finley_IndexList* index_list=NULL;
71 Finley_resetError();
72 Finley_NodeDistribution *row_degreeOfFreedomDistribution;
73 Finley_NodeDistribution *col_degreeOfFreedomDistribution;
74
75 time0=Finley_timer();
76
77 #ifdef PASO_MPI
78 if (reduce_col_order) {
79 n=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
80 colLabel=mesh->Nodes->reducedDegreeOfFreedom;
81 } else {
82 n=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
83 colLabel=mesh->Nodes->degreeOfFreedom;
84 }
85
86 if (reduce_row_order) {
87 n=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
88 rowLabel=mesh->Nodes->reducedDegreeOfFreedom;
89 } else {
90 n=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
91 rowLabel=mesh->Nodes->degreeOfFreedom;
92 }
93 #else
94 if (reduce_col_order) {
95 n=mesh->Nodes->reducedNumDegreesOfFreedom;
96 colLabel=mesh->Nodes->reducedDegreeOfFreedom;
97 } else {
98 n=mesh->Nodes->numDegreesOfFreedom;
99 colLabel=mesh->Nodes->degreeOfFreedom;
100 }
101
102 if (reduce_row_order) {
103 n=mesh->Nodes->reducedNumDegreesOfFreedom;
104 rowLabel=mesh->Nodes->reducedDegreeOfFreedom;
105 } else {
106 n=mesh->Nodes->numDegreesOfFreedom;
107 rowLabel=mesh->Nodes->degreeOfFreedom;
108 }
109 #endif
110
111 index_list=TMPMEMALLOC(n,Finley_IndexList);
112 ptr=MEMALLOC(n+1,index_t);
113 if (! (Finley_checkPtr(index_list) || Finley_checkPtr(ptr)) ) {
114 #pragma omp parallel private(i,s,itmp)
115 {
116 #pragma omp for schedule(static)
117 for(i=0;i<n;++i) {
118 index_list[i].extension=NULL;
119 index_list[i].n=0;
120 }
121 /* insert contributions from element matrices into colums index index_list: */
122 printf("ksteube mesh->Elements\n");
123 Finley_IndexList_insertElements(index_list,mesh,mesh->Elements,
124 reduce_row_order,rowLabel,reduce_col_order,colLabel);
125 printf("ksteube mesh->FaceElements\n");
126 Finley_IndexList_insertElements(index_list,mesh,mesh->FaceElements,
127 reduce_row_order,rowLabel,reduce_col_order,colLabel);
128 printf("ksteube mesh->ContactElements\n");
129 Finley_IndexList_insertElements(index_list,mesh,mesh->ContactElements,
130 reduce_row_order,rowLabel,reduce_col_order,colLabel);
131 printf("ksteube mesh->Points\n");
132 Finley_IndexList_insertElements(index_list,mesh,mesh->Points,
133 reduce_row_order,rowLabel,reduce_col_order,colLabel);
134 printf("ksteube done with 4 calls to Finley_IndexList_insertElements\n");
135 /* get the number of connections per row */
136 #pragma omp for schedule(static)
137 for(i=0;i<n;++i) {
138 ptr[i]=Finley_IndexList_count(&index_list[i]);
139 }
140 /* accumalate ptr */
141 /* OMP */
142 #pragma omp master
143 {
144 s=0;
145 for(i=0;i<n;++i) {
146 itmp=ptr[i];
147 ptr[i]=s;
148 s+=itmp;
149 }
150 ptr[n]=s;
151 /* fill index */
152 index=MEMALLOC(ptr[n],index_t);
153 }
154 #pragma omp barrier
155 if (! Finley_checkPtr(index)) {
156 #pragma omp for schedule(static)
157 for(i=0;i<n;++i) Finley_IndexList_toArray(&index_list[i],&index[ptr[i]]);
158 }
159 }
160 }
161 /* all done */
162 /* clean up */
163 if (index_list!=NULL) {
164 #pragma omp parallel for private(i)
165 for(i=0;i<n;++i) Finley_IndexList_free(index_list[i].extension);
166 }
167 TMPMEMFREE(index_list);
168 #ifdef Finley_TRACE
169 printf("timing: mesh to matrix pattern: %.4e sec\n",Finley_timer()-time0);
170 #endif
171 if (! Finley_noError()) {
172 MEMFREE(ptr);
173 MEMFREE(index);
174 return NULL;
175 } else {
176 Paso_SystemMatrixPattern *pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,n,ptr,index,mesh->MPIInfo); /* add distr info */
177 // Use a getReference method to get the DOF distributions to avoid memory leaks
178 if (reduce_col_order) {
179 col_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
180 }
181 else {
182 col_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
183 }
184 if (reduce_row_order) {
185 row_degreeOfFreedomDistribution = mesh->Nodes->reducedDegreeOfFreedomDistribution;
186 }
187 else {
188 row_degreeOfFreedomDistribution = mesh->Nodes->degreeOfFreedomDistribution;
189 }
190 pattern->row_degreeOfFreedomDistribution = row_degreeOfFreedomDistribution;
191 pattern->col_degreeOfFreedomDistribution = col_degreeOfFreedomDistribution;
192 return pattern;
193 }
194 }
195 /*
196 * $Log$
197 * Revision 1.5 2005/09/15 03:44:22 jgs
198 * Merge of development branch dev-02 back to main trunk on 2005-09-15
199 *
200 * Revision 1.4 2005/09/01 03:31:35 jgs
201 * Merge of development branch dev-02 back to main trunk on 2005-09-01
202 *
203 * Revision 1.3.2.2 2005/09/07 06:26:19 gross
204 * the solver from finley are put into the standalone package paso now
205 *
206 * Revision 1.3.2.1 2005/08/24 02:02:18 gross
207 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
208 *
209 * Revision 1.3 2005/07/08 04:07:52 jgs
210 * Merge of development branch back to main trunk on 2005-07-08
211 *
212 * Revision 1.2 2004/12/15 07:08:33 jgs
213 * *** empty log message ***
214 *
215 * Revision 1.1.2.5 2005/06/30 01:53:55 gross
216 * a bug in coloring fixed
217 *
218 * Revision 1.1.2.4 2005/06/29 02:34:51 gross
219 * some changes towards 64 integers in finley
220 *
221 * Revision 1.1.2.3 2004/11/24 01:37:14 gross
222 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
223 *
224 * Revision 1.1.2.1 2004/11/15 01:01:09 gross
225 * and anotther missing file
226 *
227 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26