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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 415 - (show annotations)
Wed Jan 4 05:37:33 2006 UTC (14 years, 4 months ago) by gross
File MIME type: text/plain
File size: 6431 byte(s)
a better way of representing the matrix format type is introduced. this is needed for the Paradiso and UMFPACK interface
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 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 if ((type & MATRIX_FORMAT_BLK1) && (row_block_size!=col_block_size || col_block_size>3) ) {
47 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 out->row_block_size=1;
53 out->col_block_size=1;
54 } 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 out->row_block_size=row_block_size;
61 out->col_block_size=col_block_size;
62 }
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 if ((type & MATRIX_FORMAT_BLK1) && (row_block_size!=col_block_size || col_block_size>3) ) {
73 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 }
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