paso/src/SystemMatrix.c

/* $Id$ */

/**************************************************************/

/* Paso: SystemMatrix */

/**************************************************************/

/* Copyrights by ACcESS Australia 2003, 2004,2005 */
/* Author: gross@access.edu.au */

/**************************************************************/

#include "Paso.h"
#include "SystemMatrix.h"

/**************************************************************/

/* allocates a SystemMatrix of type type using the given matrix pattern 
   if type is UNKOWN CSR is used.
   if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
   if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
   Values are initialized by zero.  */

Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size) {
  double time0;
  Paso_SystemMatrix*out=NULL;
  Paso_resetError();
  time0=Paso_timer();
  dim_t n_norm,i;
  out=MEMALLOC(1,Paso_SystemMatrix);
  if (! Paso_checkPtr(out)) {  
     out->pattern=NULL;  
     out->direct=NULL;  
     out->iterative=NULL;
     out->val=NULL;  
     out->reference_counter=1;
     /* check the matrix type */
     switch(type) {
        case CSC:
          out->type=CSC;
          if (row_block_size!=col_block_size || col_block_size>3) {
             out->row_block_size=1;
             out->col_block_size=1;
             out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,col_block_size,row_block_size);
          } else {
             out->pattern=Paso_SystemMatrixPattern_reference(pattern);
             out->row_block_size=row_block_size;
             out->col_block_size=col_block_size;
          }
          break;
        case CSR:
           out->type=CSR;
           if (row_block_size!=col_block_size || col_block_size>3) {
              out->row_block_size=1;
              out->col_block_size=1;
              out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,row_block_size,col_block_size);
           } else { 
              out->pattern=Paso_SystemMatrixPattern_reference(pattern);
              out->row_block_size=row_block_size;
              out->col_block_size=col_block_size;
          }
          break;
        case CSC_BLK1:
          out->type=CSC;
          out->row_block_size=1;
          out->col_block_size=1;
          if (row_block_size==1 && col_block_size==1) {
              out->pattern=Paso_SystemMatrixPattern_reference(pattern);
          } else {
             out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,col_block_size,row_block_size);
          }
          break;
        case CSR_BLK1:
          out->type=CSR;
          out->row_block_size=1;
          out->col_block_size=1;
          if (row_block_size==1 && col_block_size==1) {
              out->pattern=Paso_SystemMatrixPattern_reference(pattern);
          } else {
             out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,row_block_size,col_block_size);
          }
          break;
        case CSC_SYM:
        case CSC_BLK1_SYM:
          out->type=CSC_SYM;
          Paso_setError(TYPE_ERROR,"convertion of matrix pattern for symmetric CSC is not implemented yet.");
          return NULL;
        default:
          Paso_setError(TYPE_ERROR,"unknown matrix type identifier.");
          return NULL;
     }
     if (out->type==CSC || out->type==CSC_SYM ) {
         out->num_rows=out->pattern->n_index;
         out->num_cols=out->pattern->n_ptr;
         n_norm = out->num_cols * out->col_block_size;

     } else {
         out->num_rows=out->pattern->n_ptr;
         out->num_cols=out->pattern->n_index;
         n_norm = out->num_rows * out->row_block_size;
     } 
     out->logical_row_block_size=row_block_size;
     out->logical_col_block_size=col_block_size;
     out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
     out->block_size=out->row_block_size*out->col_block_size;
     out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
     /* allocate memory for matrix entries */
     out->val=MEMALLOC(out->len,double);
     out->normalizer=MEMALLOC(n_norm,double);
     out->normalizer_is_valid=FALSE;
     if (! Paso_checkPtr(out->val)) {
        Paso_SystemMatrix_setValues(out,DBLE(0));
     }
     if (! Paso_checkPtr(out->normalizer)) {
         #pragma omp parallel for private(i) schedule(static)
         for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
     }
  }
  /* all done: */
  if (! Paso_noError()) {
    Paso_SystemMatrix_dealloc(out);
    return NULL;
  } else {
    #ifdef Paso_TRACE
    printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
    printf("Paso_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);
    #endif
    return out;
  }
}

/* returns a reference to Paso_SystemMatrix in */

Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
   if (in!=NULL) ++(in->reference_counter);
   return NULL;
}

/* deallocates a SystemMatrix: */

void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {
  if (in!=NULL) {
     in->reference_counter--;
     if (in->reference_counter<=0) {
        MEMFREE(in->val);
        MEMFREE(in->normalizer);
        Paso_SystemMatrixPattern_dealloc(in->pattern);
        Paso_solve_free(in); 
        MEMFREE(in);
        #ifdef Paso_TRACE
        printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");
        #endif
     }
   }
}
/*
 * $Log$
 * Revision 1.2  2005/09/15 03:44:38  jgs
 * Merge of development branch dev-02 back to main trunk on 2005-09-15
 *
 * Revision 1.1.2.2  2005/09/07 00:59:08  gross
 * some inconsistent renaming fixed to make the linking work.
 *
 * Revision 1.1.2.1  2005/09/05 06:29:47  gross
 * These files have been extracted from finley to define a stand alone libray for iterative
 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
 * has not been tested yet.
 *
 *
 */
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	*/
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision