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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1028 - (hide annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 7084 byte(s)
modifications to be compliant with _WIN32. The substitutes for asinh, acosh, atanh are still missing (erf will through an exception)
1 jgs 150 /* $Id$ */
2    
3 dhawcroft 631
4     /*
5     ********************************************************************************
6 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
7 dhawcroft 631 * *
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 jgs 150 /**************************************************************/
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 gross 432
39 jgs 150 double time0;
40     Paso_SystemMatrix*out=NULL;
41 gross 1028 dim_t n_norm,i;
42    
43 jgs 150 Paso_resetError();
44     time0=Paso_timer();
45     out=MEMALLOC(1,Paso_SystemMatrix);
46     if (! Paso_checkPtr(out)) {
47     out->pattern=NULL;
48 gross 425 out->solver_package=PASO_PASO;
49     out->solver=NULL;
50 jgs 150 out->val=NULL;
51     out->reference_counter=1;
52 gross 415 out->type=type;
53     if (type & MATRIX_FORMAT_CSC) {
54     if (type & MATRIX_FORMAT_SYM) {
55     Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
56     return NULL;
57     } else {
58 gross 425 if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
59 gross 415 if (type & MATRIX_FORMAT_OFFSET1) {
60     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,col_block_size,row_block_size);
61     } else {
62     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,col_block_size,row_block_size);
63     }
64 jgs 150 out->row_block_size=1;
65     out->col_block_size=1;
66 gross 415 } else {
67     if ( (type & MATRIX_FORMAT_OFFSET1) ==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
68     out->pattern=Paso_SystemMatrixPattern_reference(pattern);
69     } else {
70     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT,1,1);
71     }
72 jgs 150 out->row_block_size=row_block_size;
73     out->col_block_size=col_block_size;
74 gross 415 }
75     }
76     out->num_rows=out->pattern->n_index;
77     out->num_cols=out->pattern->n_ptr;
78     n_norm = out->num_cols * out->col_block_size;
79     } else {
80     if (type & MATRIX_FORMAT_SYM) {
81     Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
82     return NULL;
83     } else {
84 gross 425 if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
85 gross 415 if (type & MATRIX_FORMAT_OFFSET1) {
86     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,row_block_size,col_block_size);
87     } else {
88     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,row_block_size,col_block_size);
89     }
90     out->row_block_size=1;
91     out->col_block_size=1;
92     } else {
93     if ((type & MATRIX_FORMAT_OFFSET1)==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
94     out->pattern=Paso_SystemMatrixPattern_reference(pattern);
95     } else {
96     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT,1,1);
97     }
98     out->row_block_size=row_block_size;
99     out->col_block_size=col_block_size;
100     }
101     }
102     out->num_rows=out->pattern->n_index;
103     out->num_cols=out->pattern->n_ptr;
104     n_norm = out->num_cols * out->col_block_size;
105     out->num_rows=out->pattern->n_ptr;
106     out->num_cols=out->pattern->n_index;
107     n_norm = out->num_rows * out->row_block_size;
108 jgs 150 }
109     out->logical_row_block_size=row_block_size;
110     out->logical_col_block_size=col_block_size;
111     out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
112     out->block_size=out->row_block_size*out->col_block_size;
113     out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
114     /* allocate memory for matrix entries */
115     out->val=MEMALLOC(out->len,double);
116     out->normalizer=MEMALLOC(n_norm,double);
117     out->normalizer_is_valid=FALSE;
118     if (! Paso_checkPtr(out->val)) {
119     Paso_SystemMatrix_setValues(out,DBLE(0));
120     }
121     if (! Paso_checkPtr(out->normalizer)) {
122     #pragma omp parallel for private(i) schedule(static)
123     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
124     }
125     }
126     /* all done: */
127     if (! Paso_noError()) {
128     Paso_SystemMatrix_dealloc(out);
129     return NULL;
130     } else {
131     #ifdef Paso_TRACE
132     printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
133     printf("Paso_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);
134     #endif
135     return out;
136     }
137     }
138    
139     /* returns a reference to Paso_SystemMatrix in */
140    
141     Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
142     if (in!=NULL) ++(in->reference_counter);
143     return NULL;
144     }
145    
146     /* deallocates a SystemMatrix: */
147    
148     void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {
149     if (in!=NULL) {
150     in->reference_counter--;
151     if (in->reference_counter<=0) {
152     MEMFREE(in->val);
153     MEMFREE(in->normalizer);
154     Paso_SystemMatrixPattern_dealloc(in->pattern);
155     Paso_solve_free(in);
156     MEMFREE(in);
157     #ifdef Paso_TRACE
158     printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");
159     #endif
160     }
161     }
162     }
163     /*
164     * $Log$
165     * Revision 1.2 2005/09/15 03:44:38 jgs
166     * Merge of development branch dev-02 back to main trunk on 2005-09-15
167     *
168     * Revision 1.1.2.2 2005/09/07 00:59:08 gross
169     * some inconsistent renaming fixed to make the linking work.
170     *
171     * Revision 1.1.2.1 2005/09/05 06:29:47 gross
172     * These files have been extracted from finley to define a stand alone libray for iterative
173     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
174     * has not been tested yet.
175     *
176     *
177     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26