/[escript]/trunk/paso/src/SystemMatrix.c
ViewVC logotype

Annotation of /trunk/paso/src/SystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1562 - (hide annotations)
Wed May 21 13:04:40 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/plain
File size: 8717 byte(s)
The algebraic upwinding with MPI. The case of boundary constraint needs still some attention. 


1 ksteube 1312
2 jgs 150 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 dhawcroft 631
16 jgs 150 /**************************************************************/
17    
18     /* Paso: SystemMatrix */
19    
20     /**************************************************************/
21    
22     /* Author: gross@access.edu.au */
23    
24     /**************************************************************/
25    
26     #include "SystemMatrix.h"
27    
28     /**************************************************************/
29    
30     /* allocates a SystemMatrix of type type using the given matrix pattern
31     if type is UNKOWN CSR is used.
32     if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
33     if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
34     Values are initialized by zero. */
35    
36     Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size) {
37 gross 432
38 jgs 150 double time0;
39     Paso_SystemMatrix*out=NULL;
40 gross 1028 dim_t n_norm,i;
41 ksteube 1312 Paso_SystemMatrixType pattern_format_out;
42 gross 1028
43 jgs 150 Paso_resetError();
44     time0=Paso_timer();
45     out=MEMALLOC(1,Paso_SystemMatrix);
46     if (! Paso_checkPtr(out)) {
47 ksteube 1312 out->type=type;
48 jgs 150 out->pattern=NULL;
49 ksteube 1312 out->row_distribution=NULL;
50     out->col_distribution=NULL;
51     out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
52 gross 1552 out->row_coupler=NULL;
53     out->col_coupler=NULL;
54 ksteube 1312 out->mainBlock=NULL;
55 gross 1552 out->row_coupleBlock=NULL;
56     out->col_coupleBlock=NULL;
57 ksteube 1312 out->normalizer_is_valid=FALSE;
58     out->normalizer=NULL;
59 gross 425 out->solver_package=PASO_PASO;
60     out->solver=NULL;
61 ksteube 1312 out->trilinos_data=NULL;
62 jgs 150 out->reference_counter=1;
63 ksteube 1312
64 gross 1552
65 ksteube 1312 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
66     /* ====== compressed sparse columns === */
67 gross 415 if (type & MATRIX_FORMAT_CSC) {
68     if (type & MATRIX_FORMAT_SYM) {
69     Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
70     } else {
71 gross 425 if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
72 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
73 jgs 150 out->row_block_size=1;
74     out->col_block_size=1;
75 gross 415 } else {
76 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
77 jgs 150 out->row_block_size=row_block_size;
78     out->col_block_size=col_block_size;
79 gross 415 }
80     }
81 ksteube 1312 out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
82     out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
83 gross 415 } else {
84 ksteube 1312 /* ====== compressed sparse row === */
85 gross 415 if (type & MATRIX_FORMAT_SYM) {
86     Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
87     } else {
88 gross 425 if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
89 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
90 gross 415 out->row_block_size=1;
91     out->col_block_size=1;
92     } else {
93 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
94 gross 415 out->row_block_size=row_block_size;
95     out->col_block_size=col_block_size;
96     }
97     }
98 ksteube 1312 out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
99     out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
100 jgs 150 }
101     out->logical_row_block_size=row_block_size;
102     out->logical_col_block_size=col_block_size;
103     out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
104     out->block_size=out->row_block_size*out->col_block_size;
105 gross 1552 out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);
106     out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);
107 ksteube 1312 /* this should be bypassed if trilinos is used */
108     if (type & MATRIX_FORMAT_TRILINOS_CRS) {
109     #ifdef TRILINOS
110     out->trilinos_data=Paso_TRILINOS_alloc();
111     #endif
112     } else {
113     out->solver_package=PASO_PASO;
114     out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size);
115 gross 1552 out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size);
116     out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size);
117 ksteube 1312 /* allocate memory for matrix entries */
118     if (type & MATRIX_FORMAT_CSC) {
119     n_norm = out->mainBlock->numCols * out->col_block_size;
120     } else {
121     n_norm = out->mainBlock->numRows * out->row_block_size;
122     }
123     out->normalizer=MEMALLOC(n_norm,double);
124     out->normalizer_is_valid=FALSE;
125     if (! Paso_checkPtr(out->normalizer)) {
126     #pragma omp parallel for private(i) schedule(static)
127     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
128     }
129     }
130 jgs 150 }
131     /* all done: */
132     if (! Paso_noError()) {
133 ksteube 1312 Paso_SystemMatrix_free(out);
134 jgs 150 return NULL;
135     } else {
136     #ifdef Paso_TRACE
137     printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
138 ksteube 1312 printf("Paso_SystemMatrix_alloc: system matrix has been allocated.\n");
139 jgs 150 #endif
140     return out;
141     }
142     }
143    
144     /* returns a reference to Paso_SystemMatrix in */
145    
146     Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
147     if (in!=NULL) ++(in->reference_counter);
148     return NULL;
149     }
150    
151     /* deallocates a SystemMatrix: */
152    
153 ksteube 1312 void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
154 jgs 150 if (in!=NULL) {
155     in->reference_counter--;
156     if (in->reference_counter<=0) {
157 ksteube 1312 Paso_SystemMatrixPattern_free(in->pattern);
158     Paso_Distribution_free(in->row_distribution);
159     Paso_Distribution_free(in->col_distribution);
160     Paso_MPIInfo_free(in->mpi_info);
161 gross 1552 Paso_Coupler_free(in->row_coupler);
162     Paso_Coupler_free(in->col_coupler);
163 ksteube 1312 Paso_SparseMatrix_free(in->mainBlock);
164 gross 1552 Paso_SparseMatrix_free(in->col_coupleBlock);
165     Paso_SparseMatrix_free(in->row_coupleBlock);
166 jgs 150 MEMFREE(in->normalizer);
167     Paso_solve_free(in);
168 ksteube 1312 #ifdef TRILINOS
169     Paso_TRILINOS_free(in->trilinos_data);
170     #endif
171 jgs 150 MEMFREE(in);
172     #ifdef Paso_TRACE
173 ksteube 1312 printf("Paso_SystemMatrix_free: system matrix as been deallocated.\n");
174 jgs 150 #endif
175     }
176     }
177     }
178 gross 1562 void Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,double* in)
179 ksteube 1312 {
180 gross 1552 Paso_SystemMatrix_startColCollect(A,in);
181 ksteube 1312 }
182     double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
183     {
184 gross 1552 return Paso_SystemMatrix_finishColCollect(A);
185 ksteube 1312 }
186    
187 gross 1562 void Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,double* in)
188 gross 1552 {
189     Paso_Coupler_startCollect(A->col_coupler, in);
190     }
191     double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A)
192     {
193     Paso_Coupler_finishCollect(A->col_coupler);
194     return A->col_coupler->recv_buffer;
195     }
196 gross 1562 void Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,double* in)
197 gross 1552 {
198     Paso_Coupler_startCollect(A->row_coupler, in);
199     }
200     double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A)
201     {
202     Paso_Coupler_finishCollect(A->row_coupler);
203     return A->row_coupler->recv_buffer;
204     }
205    
206 gross 1407 dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
207 ksteube 1312 return A->mainBlock->numRows * A->row_block_size;
208     }
209    
210 gross 1407 dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
211 ksteube 1312 return A->mainBlock->numCols * A->col_block_size;
212     }
213     dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {
214     if (A->type & MATRIX_FORMAT_CSC) {
215     return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
216     } else {
217     return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
218     }
219     }
220     dim_t Paso_SystemMatrix_getGlobalNumCols(Paso_SystemMatrix* A) {
221     if (A->type & MATRIX_FORMAT_CSC) {
222     return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
223     } else {
224     return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
225     }
226    
227     }
228 gross 1361 dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
229     return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
230     }
231 ksteube 1312

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26