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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 969 - (hide annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 9 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 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 102
13     /**************************************************************/
14    
15     /* Finley: Mesh */
16    
17     /**************************************************************/
18    
19 jgs 150 /* Author: gross@access.edu.au */
20     /* Version: $Id$ */
21 jgs 102
22     /**************************************************************/
23    
24 jgs 150 #include "Mesh.h"
25 jgs 102 #include "IndexList.h"
26    
27     /**************************************************************/
28    
29     /* returns a reference to the matrix pattern */
30    
31 jgs 150 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 jgs 102 /* 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 jgs 150 if (Finley_noError()) {
49 jgs 102 if (reduce_row_order) {
50     if (reduce_col_order) {
51 jgs 150 out=Paso_SystemMatrixPattern_reference(mesh->ReducedReducedPattern);
52 jgs 102 } else {
53 jgs 150 out=Paso_SystemMatrixPattern_reference(mesh->ReducedFullPattern);
54 jgs 102 }
55     } else {
56     if (reduce_col_order) {
57 jgs 150 out=Paso_SystemMatrixPattern_reference(mesh->FullReducedPattern);
58 jgs 102 } else {
59 jgs 150 out=Paso_SystemMatrixPattern_reference(mesh->FullFullPattern);
60 jgs 102 }
61     }
62     }
63 ksteube 969 out->MPIInfo = Paso_MPIInfo_getReference( mesh->MPIInfo );
64 jgs 102 return out;
65     }
66 jgs 150 Paso_SystemMatrixPattern* Finley_makePattern(Finley_Mesh *mesh,bool_t reduce_row_order, bool_t reduce_col_order) {
67 jgs 102 double time0;
68 jgs 123 dim_t i,n;
69     index_t s,itmp,*rowLabel=NULL,*colLabel=NULL, *ptr=NULL, *index=NULL;
70 jgs 102 Finley_IndexList* index_list=NULL;
71 jgs 150 Finley_resetError();
72 ksteube 969 Finley_NodeDistribution *row_degreeOfFreedomDistribution;
73     Finley_NodeDistribution *col_degreeOfFreedomDistribution;
74    
75 jgs 102 time0=Finley_timer();
76 ksteube 969
77     #ifdef PASO_MPI
78 jgs 102 if (reduce_col_order) {
79 ksteube 969 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 jgs 102 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 ksteube 969 #endif
110 jgs 102
111     index_list=TMPMEMALLOC(n,Finley_IndexList);
112 jgs 123 ptr=MEMALLOC(n+1,index_t);
113 jgs 102 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 ksteube 969 printf("ksteube mesh->Elements\n");
123     Finley_IndexList_insertElements(index_list,mesh,mesh->Elements,
124 jgs 102 reduce_row_order,rowLabel,reduce_col_order,colLabel);
125 ksteube 969 printf("ksteube mesh->FaceElements\n");
126     Finley_IndexList_insertElements(index_list,mesh,mesh->FaceElements,
127 jgs 102 reduce_row_order,rowLabel,reduce_col_order,colLabel);
128 ksteube 969 printf("ksteube mesh->ContactElements\n");
129     Finley_IndexList_insertElements(index_list,mesh,mesh->ContactElements,
130 jgs 102 reduce_row_order,rowLabel,reduce_col_order,colLabel);
131 ksteube 969 printf("ksteube mesh->Points\n");
132     Finley_IndexList_insertElements(index_list,mesh,mesh->Points,
133 jgs 102 reduce_row_order,rowLabel,reduce_col_order,colLabel);
134 ksteube 969 printf("ksteube done with 4 calls to Finley_IndexList_insertElements\n");
135 jgs 102 /* get the number of connections per row */
136     #pragma omp for schedule(static)
137 jgs 123 for(i=0;i<n;++i) {
138     ptr[i]=Finley_IndexList_count(&index_list[i]);
139     }
140 jgs 102 /* 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 jgs 123 index=MEMALLOC(ptr[n],index_t);
153 jgs 102 }
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 jgs 149 #ifdef Finley_TRACE
169 jgs 102 printf("timing: mesh to matrix pattern: %.4e sec\n",Finley_timer()-time0);
170 jgs 149 #endif
171 jgs 150 if (! Finley_noError()) {
172 jgs 102 MEMFREE(ptr);
173     MEMFREE(index);
174     return NULL;
175     } else {
176 ksteube 969 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 jgs 102 }
194     }
195     /*
196     * $Log$
197 jgs 150 * 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 jgs 149 * 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 jgs 150 * 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 jgs 149 * 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 jgs 123 * 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 jgs 102 * Revision 1.2 2004/12/15 07:08:33 jgs
213     * *** empty log message ***
214     *
215 jgs 123 * 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 jgs 102 * 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