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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 425 - (hide annotations)
Tue Jan 10 04:10:39 2006 UTC (13 years, 9 months ago) by gross
File MIME type: text/plain
File size: 6438 byte(s)
The sparse solver can be called by paso now. 

the building has been change to reduce some code redundancy:
now all scons SCscripts are importing scons/esys_options.py which
imports platform specific settings. 



1 jgs 150 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Paso: SystemMatrix */
6    
7     /**************************************************************/
8    
9     /* Copyrights by ACcESS Australia 2003, 2004,2005 */
10     /* Author: gross@access.edu.au */
11    
12     /**************************************************************/
13    
14     #include "Paso.h"
15     #include "SystemMatrix.h"
16    
17     /**************************************************************/
18    
19     /* allocates a SystemMatrix of type type using the given matrix pattern
20     if type is UNKOWN CSR is used.
21     if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
22     if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
23     Values are initialized by zero. */
24    
25     Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size) {
26     double time0;
27     Paso_SystemMatrix*out=NULL;
28     Paso_resetError();
29     time0=Paso_timer();
30     dim_t n_norm,i;
31     out=MEMALLOC(1,Paso_SystemMatrix);
32     if (! Paso_checkPtr(out)) {
33     out->pattern=NULL;
34 gross 425 out->solver_package=PASO_PASO;
35     out->solver=NULL;
36 jgs 150 out->val=NULL;
37     out->reference_counter=1;
38 gross 415 out->type=type;
39    
40    
41     if (type & MATRIX_FORMAT_CSC) {
42     if (type & MATRIX_FORMAT_SYM) {
43     Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
44     return NULL;
45     } else {
46 gross 425 if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
47 gross 415 if (type & MATRIX_FORMAT_OFFSET1) {
48     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,col_block_size,row_block_size);
49     } else {
50     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,col_block_size,row_block_size);
51     }
52 jgs 150 out->row_block_size=1;
53     out->col_block_size=1;
54 gross 415 } else {
55     if ( (type & MATRIX_FORMAT_OFFSET1) ==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
56     out->pattern=Paso_SystemMatrixPattern_reference(pattern);
57     } else {
58     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT,1,1);
59     }
60 jgs 150 out->row_block_size=row_block_size;
61     out->col_block_size=col_block_size;
62 gross 415 }
63     }
64     out->num_rows=out->pattern->n_index;
65     out->num_cols=out->pattern->n_ptr;
66     n_norm = out->num_cols * out->col_block_size;
67     } else {
68     if (type & MATRIX_FORMAT_SYM) {
69     Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
70     return NULL;
71     } else {
72 gross 425 if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
73 gross 415 if (type & MATRIX_FORMAT_OFFSET1) {
74     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,row_block_size,col_block_size);
75     } else {
76     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,row_block_size,col_block_size);
77     }
78     out->row_block_size=1;
79     out->col_block_size=1;
80     } else {
81     if ((type & MATRIX_FORMAT_OFFSET1)==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
82     out->pattern=Paso_SystemMatrixPattern_reference(pattern);
83     } else {
84     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT,1,1);
85     }
86     out->row_block_size=row_block_size;
87     out->col_block_size=col_block_size;
88     }
89     }
90     out->num_rows=out->pattern->n_index;
91     out->num_cols=out->pattern->n_ptr;
92     n_norm = out->num_cols * out->col_block_size;
93     out->num_rows=out->pattern->n_ptr;
94     out->num_cols=out->pattern->n_index;
95     n_norm = out->num_rows * out->row_block_size;
96 jgs 150 }
97     out->logical_row_block_size=row_block_size;
98     out->logical_col_block_size=col_block_size;
99     out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
100     out->block_size=out->row_block_size*out->col_block_size;
101     out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
102     /* allocate memory for matrix entries */
103     out->val=MEMALLOC(out->len,double);
104     out->normalizer=MEMALLOC(n_norm,double);
105     out->normalizer_is_valid=FALSE;
106     if (! Paso_checkPtr(out->val)) {
107     Paso_SystemMatrix_setValues(out,DBLE(0));
108     }
109     if (! Paso_checkPtr(out->normalizer)) {
110     #pragma omp parallel for private(i) schedule(static)
111     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
112     }
113     }
114     /* all done: */
115     if (! Paso_noError()) {
116     Paso_SystemMatrix_dealloc(out);
117     return NULL;
118     } else {
119     #ifdef Paso_TRACE
120     printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
121     printf("Paso_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);
122     #endif
123     return out;
124     }
125     }
126    
127     /* returns a reference to Paso_SystemMatrix in */
128    
129     Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
130     if (in!=NULL) ++(in->reference_counter);
131     return NULL;
132     }
133    
134     /* deallocates a SystemMatrix: */
135    
136     void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {
137     if (in!=NULL) {
138     in->reference_counter--;
139     if (in->reference_counter<=0) {
140     MEMFREE(in->val);
141     MEMFREE(in->normalizer);
142     Paso_SystemMatrixPattern_dealloc(in->pattern);
143     Paso_solve_free(in);
144     MEMFREE(in);
145     #ifdef Paso_TRACE
146     printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");
147     #endif
148     }
149     }
150     }
151     /*
152     * $Log$
153     * Revision 1.2 2005/09/15 03:44:38 jgs
154     * Merge of development branch dev-02 back to main trunk on 2005-09-15
155     *
156     * Revision 1.1.2.2 2005/09/07 00:59:08 gross
157     * some inconsistent renaming fixed to make the linking work.
158     *
159     * Revision 1.1.2.1 2005/09/05 06:29:47 gross
160     * These files have been extracted from finley to define a stand alone libray for iterative
161     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
162     * has not been tested yet.
163     *
164     *
165     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26