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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 10 months ago) by jgs
File MIME type: text/plain
File size: 5999 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 /* $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 out->direct=NULL;
35 out->iterative=NULL;
36 out->val=NULL;
37 out->reference_counter=1;
38 /* check the matrix type */
39 switch(type) {
40 case CSC:
41 out->type=CSC;
42 if (row_block_size!=col_block_size || col_block_size>3) {
43 out->row_block_size=1;
44 out->col_block_size=1;
45 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,col_block_size,row_block_size);
46 } else {
47 out->pattern=Paso_SystemMatrixPattern_reference(pattern);
48 out->row_block_size=row_block_size;
49 out->col_block_size=col_block_size;
50 }
51 break;
52 case CSR:
53 out->type=CSR;
54 if (row_block_size!=col_block_size || col_block_size>3) {
55 out->row_block_size=1;
56 out->col_block_size=1;
57 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,row_block_size,col_block_size);
58 } else {
59 out->pattern=Paso_SystemMatrixPattern_reference(pattern);
60 out->row_block_size=row_block_size;
61 out->col_block_size=col_block_size;
62 }
63 break;
64 case CSC_BLK1:
65 out->type=CSC;
66 out->row_block_size=1;
67 out->col_block_size=1;
68 if (row_block_size==1 && col_block_size==1) {
69 out->pattern=Paso_SystemMatrixPattern_reference(pattern);
70 } else {
71 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,col_block_size,row_block_size);
72 }
73 break;
74 case CSR_BLK1:
75 out->type=CSR;
76 out->row_block_size=1;
77 out->col_block_size=1;
78 if (row_block_size==1 && col_block_size==1) {
79 out->pattern=Paso_SystemMatrixPattern_reference(pattern);
80 } else {
81 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,row_block_size,col_block_size);
82 }
83 break;
84 case CSC_SYM:
85 case CSC_BLK1_SYM:
86 out->type=CSC_SYM;
87 Paso_setError(TYPE_ERROR,"convertion of matrix pattern for symmetric CSC is not implemented yet.");
88 return NULL;
89 default:
90 Paso_setError(TYPE_ERROR,"unknown matrix type identifier.");
91 return NULL;
92 }
93 if (out->type==CSC || out->type==CSC_SYM ) {
94 out->num_rows=out->pattern->n_index;
95 out->num_cols=out->pattern->n_ptr;
96 n_norm = out->num_cols * out->col_block_size;
97
98 } else {
99 out->num_rows=out->pattern->n_ptr;
100 out->num_cols=out->pattern->n_index;
101 n_norm = out->num_rows * out->row_block_size;
102 }
103 out->logical_row_block_size=row_block_size;
104 out->logical_col_block_size=col_block_size;
105 out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
106 out->block_size=out->row_block_size*out->col_block_size;
107 out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
108 /* allocate memory for matrix entries */
109 out->val=MEMALLOC(out->len,double);
110 out->normalizer=MEMALLOC(n_norm,double);
111 out->normalizer_is_valid=FALSE;
112 if (! Paso_checkPtr(out->val)) {
113 Paso_SystemMatrix_setValues(out,DBLE(0));
114 }
115 if (! Paso_checkPtr(out->normalizer)) {
116 #pragma omp parallel for private(i) schedule(static)
117 for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
118 }
119 }
120 /* all done: */
121 if (! Paso_noError()) {
122 Paso_SystemMatrix_dealloc(out);
123 return NULL;
124 } else {
125 #ifdef Paso_TRACE
126 printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
127 printf("Paso_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);
128 #endif
129 return out;
130 }
131 }
132
133 /* returns a reference to Paso_SystemMatrix in */
134
135 Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
136 if (in!=NULL) ++(in->reference_counter);
137 return NULL;
138 }
139
140 /* deallocates a SystemMatrix: */
141
142 void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {
143 if (in!=NULL) {
144 in->reference_counter--;
145 if (in->reference_counter<=0) {
146 MEMFREE(in->val);
147 MEMFREE(in->normalizer);
148 Paso_SystemMatrixPattern_dealloc(in->pattern);
149 Paso_solve_free(in);
150 MEMFREE(in);
151 #ifdef Paso_TRACE
152 printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");
153 #endif
154 }
155 }
156 }
157 /*
158 * $Log$
159 * Revision 1.2 2005/09/15 03:44:38 jgs
160 * Merge of development branch dev-02 back to main trunk on 2005-09-15
161 *
162 * Revision 1.1.2.2 2005/09/07 00:59:08 gross
163 * some inconsistent renaming fixed to make the linking work.
164 *
165 * Revision 1.1.2.1 2005/09/05 06:29:47 gross
166 * These files have been extracted from finley to define a stand alone libray for iterative
167 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
168 * has not been tested yet.
169 *
170 *
171 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26