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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 971 - (show annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 7083 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26