/[escript]/trunk/esys2/finley/src/finleyC/System.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/System.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 100 - (show annotations)
Wed Dec 15 03:48:48 2004 UTC (14 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 8669 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: SystemMatrix */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003, 2004 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "IndexList.h"
16 #include "System.h"
17 #include "Mesh.h"
18
19 /**************************************************************/
20
21 /* allocates a SystemMatrix in the CSR/Harwell boeing format from a
22 Finley_Mesh. Values are initialized by zero. */
23
24 Finley_SystemMatrix* Finley_SystemMatrix_alloc(Finley_Mesh *mesh,
25 Finley_SystemMatrixType type, int symmetric,int row_block_size, int reduce_row_order,
26 int col_block_size, int reduce_col_order) {
27 double time0,time1;
28 maybelong packed_row_block_size, packed_col_block_size, unpacked_row_block_size, unpacked_col_block_size,i,block_size,j,iptr;
29 maybelong *rowLabel=NULL,*colLabel=NULL,len_index_list=0;
30 Finley_SystemMatrix*out;
31 Finley_IndexList* index_list=NULL;
32 Finley_ErrorCode=NO_ERROR;
33
34 time0=Finley_timer();
35 /* the matrix block row_block_size x col_block_size is fully
36 incorporated into the matrix pattern. this is required for using
37 the SGI scientific library. row_offset is the offset for the row
38 index typically */
39
40 packed_row_block_size=row_block_size;
41 packed_col_block_size=col_block_size;
42 unpacked_row_block_size=1;
43 unpacked_col_block_size=1;
44
45
46 /* allocate the return value */
47
48 out=(Finley_SystemMatrix*)MEMALLOC(sizeof(Finley_SystemMatrix));
49 if (Finley_checkPtr(out)) goto clean;
50
51 if (type==UNKNOWN) {
52 out->type=FINLEY_DEFAULT_MATRIX_TYPE;
53 } else {
54 out->type=type;
55 }
56
57 out->total_row_block_size=row_block_size;
58 out->total_col_block_size=col_block_size;
59
60 out->row_block_size=unpacked_row_block_size;
61 out->col_block_size=unpacked_col_block_size;
62
63 block_size=out->row_block_size*out->col_block_size;
64
65 /* ******************************** */
66 if (reduce_col_order) {
67 out->num_cols=packed_col_block_size*mesh->Nodes->reducedNumDegreesOfFreedom;
68 colLabel=mesh->Nodes->reducedDegreeOfFreedom;
69 } else {
70 out->num_cols=packed_col_block_size*mesh->Nodes->numDegreesOfFreedom;
71 colLabel=mesh->Nodes->degreeOfFreedom;
72 }
73
74 if (reduce_row_order) {
75 out->num_rows=packed_row_block_size*mesh->Nodes->reducedNumDegreesOfFreedom;
76 rowLabel=mesh->Nodes->reducedDegreeOfFreedom;
77 } else {
78 out->num_rows=packed_row_block_size*mesh->Nodes->numDegreesOfFreedom;
79 rowLabel=mesh->Nodes->degreeOfFreedom;
80 }
81 /* ******************************** */
82
83 out->symmetric=symmetric;
84 out->index=NULL;
85 out->val=NULL;
86 out->ptr=NULL;
87 out->reference_counter=0;
88 out->lenOfVal=0;
89 out->solve=NULL;
90 out->iterative=NULL;
91
92 /* initialize the (temporary) list index_list of the colums indices: */
93 switch(out->type) {
94 case CSR:
95 len_index_list=out->num_rows;
96 break;
97 case CSC:
98 len_index_list=out->num_cols;
99 break;
100 default:
101 Finley_ErrorCode=TYPE_ERROR;
102 sprintf(Finley_ErrorMsg,"Unknown matrix type.");
103 goto clean;
104 } /* switch out->type */
105
106 index_list=(Finley_IndexList*) TMPMEMALLOC(len_index_list*sizeof(Finley_IndexList));
107 if (Finley_checkPtr(index_list)) goto clean;
108
109 #pragma omp parallel for private(i) schedule(static)
110 for(i=0;i<len_index_list;i++) {
111 index_list[i].extension=NULL;
112 index_list[i].n=0;
113 }
114
115 /* insert contributions from element matrices into colums index index_list: */
116 Finley_IndexList_insertElements(index_list,mesh->Elements,
117 reduce_row_order,packed_row_block_size,rowLabel,
118 reduce_col_order,packed_col_block_size,colLabel,
119 symmetric, out->type);
120 Finley_IndexList_insertElements(index_list,mesh->FaceElements,
121 reduce_row_order,packed_row_block_size,rowLabel,
122 reduce_col_order,packed_col_block_size,colLabel,
123 symmetric, out->type);
124 Finley_IndexList_insertElements(index_list,mesh->ContactElements,
125 reduce_row_order,packed_row_block_size,rowLabel,
126 reduce_col_order,packed_col_block_size,colLabel,
127 symmetric, out->type);
128 Finley_IndexList_insertElements(index_list,mesh->Points,
129 reduce_row_order,packed_row_block_size,rowLabel,
130 reduce_col_order,packed_col_block_size,colLabel,
131 symmetric, out->type);
132 if (Finley_ErrorCode!=NO_ERROR) goto clean;
133
134 /* allocate the ptr: */
135
136 out->ptr=(maybelong*)MEMALLOC(((len_index_list)+1)*sizeof(maybelong));
137 if (Finley_checkPtr(out->ptr)) {
138 Finley_SystemMatrix_dealloc(out);
139 goto clean;
140 }
141 #pragma omp parallel for private(i) schedule(static)
142 for(i=0;i<len_index_list;i++) out->ptr[i]=0;
143 out->ptr[len_index_list]=0;
144
145 /* count entries in each column and store to pointer vector: */
146 out->ptr[0]=PTR_OFFSET;
147 /* OMP */
148 for(i=0;i<len_index_list;i++)
149 out->ptr[i+1]=out->ptr[i]+Finley_IndexList_count(&index_list[i]);
150
151 /* allocate index and val: */
152
153 out->index=(maybelong*)MEMALLOC((out->ptr[len_index_list]-PTR_OFFSET)*sizeof(maybelong));
154 out->lenOfVal=(out->ptr[len_index_list]-PTR_OFFSET)*out->row_block_size*out->col_block_size;
155 out->val=(double*)MEMALLOC(out->lenOfVal*sizeof(double));
156 if (Finley_checkPtr(out->index) || Finley_checkPtr(out->val) ) {
157 Finley_SystemMatrix_dealloc(out);
158 goto clean;
159 }
160 #pragma omp parallel firstprivate(out,index_list,len_index_list)
161 {
162
163 #pragma omp master
164 {
165 time1=Finley_timer();
166 }
167 #pragma omp for private(i,iptr,j) schedule(static)
168 for (i=0;i< len_index_list;i++) {
169 for (iptr=out->ptr[i]-PTR_OFFSET;iptr<out->ptr[i+1]-PTR_OFFSET; iptr++) {
170 out->index[iptr]=0;
171 for (j=0;j<block_size;j++) out->val[iptr*block_size+j]=0.;
172 }
173 }
174 #pragma omp master
175 {
176 printf("timing: matrix pattern: instantation: %.4e sec\n",Finley_timer()-time1);
177 }
178
179
180 /* change the list of the row indicies into an array: */
181 /* each index is ordered by increased size */
182
183 #pragma omp for private(i) schedule(static)
184 for(i=0;i<len_index_list;i++) {
185 Finley_IndexList_toArray(&index_list[i],&(out->index[out->ptr[i]-PTR_OFFSET]));
186 qsort(&(out->index[out->ptr[i]-PTR_OFFSET]),(int)(out->ptr[i+1]-out->ptr[i]),sizeof(maybelong), Finley_comparIndex);
187 }
188
189 }
190 /* out->val is initialized to zero: */
191
192 time1=Finley_timer();
193 Finley_SystemMatrix_setValues(out,DBLE(0));
194 printf("timing: matrix pattern: initialize: %.4e sec\n",Finley_timer()-time1);
195 /* all done: */
196
197
198 clean:
199 if (index_list!=NULL) {
200 #pragma omp parallel for private(i)
201 for(i=0;i<len_index_list;i++) Finley_IndexList_free(index_list[i].extension);
202 }
203 TMPMEMFREE(index_list);
204
205 printf("timing: matrix pattern: %.4e sec\n",Finley_timer()-time0);
206
207 if (Finley_ErrorCode!=NO_ERROR) {
208 return NULL;
209 } else {
210 #ifdef Finley_TRACE
211 printf("Finley_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);
212 #endif
213 out->reference_counter++;
214 return out;
215 }
216 }
217
218 /* deallocates a SystemMatrix: */
219
220 void Finley_SystemMatrix_dealloc(Finley_SystemMatrix* in) {
221 if (in!=NULL) {
222 in->reference_counter--;
223 if (in->reference_counter<=0) {
224 MEMFREE(in->ptr);
225 MEMFREE(in->index);
226 MEMFREE(in->val);
227 Finley_SystemMatrix_solve_free(in);
228 Finley_SystemMatrix_iterative_free(in);
229 MEMFREE(in);
230 #ifdef Finley_TRACE
231 printf("Finley_SystemMatrix_dealloc: system matrix as been deallocated.\n");
232 #endif
233 }
234 }
235 }
236 /* *************************************************************/
237
238 /* some routines which help to get the matrix pattern from elements: */
239
240 /* this routine is used by qsort called in Finley_SystemMatrix_alloc */
241
242 int Finley_comparIndex(const void *index1,const void *index2){
243 maybelong Iindex1,Iindex2;
244 Iindex1=*(maybelong*)index1;
245 Iindex2=*(maybelong*)index2;
246 if (Iindex1<Iindex2) {
247 return -1;
248 } else {
249 if (Iindex1>Iindex2) {
250 return 1;
251 } else {
252 return 0;
253 }
254 }
255 }
256 /*
257 * $Log$
258 * Revision 1.3 2004/12/15 03:48:46 jgs
259 * *** empty log message ***
260 *
261 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
262 * initial import of project esys2
263 *
264 * Revision 1.1 2004/07/02 04:21:13 gross
265 * Finley C code has been included
266 *
267 *
268 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26