paso/src/SystemMatrix.c

/* $Id$ */


/*
********************************************************************************
*               Copyright   2006 by ACcESS MNRF                                *
*                                                                              * 
*                 http://www.access.edu.au                                     *
*           Primary Business: Queensland, Australia                            *
*     Licensed under the Open Software License version 3.0             *
*        http://www.opensource.org/licenses/osl-3.0.php                        *
********************************************************************************
*/

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

/* Paso: SystemMatrix */

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

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

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

#include "Paso.h"
#include "Paso_MPI.h"
#include "SystemMatrix.h"
#include "TRILINOS.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->solver_package=PASO_PASO;  
     out->solver=NULL;  
     out->val=NULL;  
     out->trilinos_data=NULL;
     out->reference_counter=1;
     out->type=type;
     /* ====== compressed sparse columns === */
     if (type & MATRIX_FORMAT_CSC) {
        if (type & MATRIX_FORMAT_SYM) {
           Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
           return NULL;
        } else {
           if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
              if (type & MATRIX_FORMAT_OFFSET1) {
                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,col_block_size,row_block_size);
              } else {
                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,col_block_size,row_block_size);
              }
              out->row_block_size=1;
              out->col_block_size=1;
           } else {
              if ( (type & MATRIX_FORMAT_OFFSET1) ==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
                  out->pattern=Paso_SystemMatrixPattern_reference(pattern);
              } else {
                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT,1,1);
              }
              out->row_block_size=row_block_size;
              out->col_block_size=col_block_size;
           }
        }
        out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
        out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
     } else {
     /* ====== compressed sparse row === */
        if (type & MATRIX_FORMAT_SYM) {
           Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
           return NULL;
        } else {
           if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3)  {
              if (type & MATRIX_FORMAT_OFFSET1) {
                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,row_block_size,col_block_size);
              } else {
                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,row_block_size,col_block_size);
              }
              out->row_block_size=1;
              out->col_block_size=1;
           } else {
              if ((type & MATRIX_FORMAT_OFFSET1)==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
                  out->pattern=Paso_SystemMatrixPattern_reference(pattern);
              } else {
                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT,1,1);
              }
              out->row_block_size=row_block_size;
              out->col_block_size=col_block_size;
           }
        }
        out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
        out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
     }
     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->myLen=(size_t)(out->pattern->myLen)*(size_t)(out->block_size);

     out->numRows = out->row_distribution->myNumComponents;
     out->myNumRows = out->row_distribution->numComponents;
     out->numCols = out->col_distribution->myNumComponents;
     out->myNumCols = out->col_distribution->numComponents;
     out->mpi_info = Paso_MPIInfo_getReference(out->pattern->mpi_info);
     /* allocate memory for matrix entries */
     if (type & MATRIX_FORMAT_TRILINOS_CRS) {
         Paso_TRILINOS_alloc(out->trilinos_data, out->pattern,out->row_block_size,out->col_block_size);
     } else {
         if (type & MATRIX_FORMAT_CSC) {
            n_norm = out->myNumCols * out->col_block_size;
         } else {
            n_norm = out->myNumRows * out->row_block_size;
         }
         out->val=MEMALLOC(out->myLen,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->numRows,(long)out->numCols);
    #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) {
        Paso_SystemMatrixPattern_dealloc(in->pattern);
        Paso_Distribution_free(in->row_distribution);
        Paso_Distribution_free(in->col_distribution);
        Paso_MPIInfo_dealloc(in->mpi_info);
        MEMFREE(in->val);
        MEMFREE(in->normalizer);
        Paso_TRILINOS_free(in->trilinos_data);
        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	********************************************************************************
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 "Paso_MPI.h"
28	#include "SystemMatrix.h"
29	#include "TRILINOS.h"
30
31	/**************************************************************/
32
33	/* allocates a SystemMatrix of type type using the given matrix pattern
34	if type is UNKOWN CSR is used.
35	if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
36	if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
37	Values are initialized by zero. */
38
39	Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size) {
40
41	double time0;
42	Paso_SystemMatrix*out=NULL;
43	Paso_resetError();
44	time0=Paso_timer();
45	dim_t n_norm,i;
46	out=MEMALLOC(1,Paso_SystemMatrix);
47	if (! Paso_checkPtr(out)) {
48	out->pattern=NULL;
49	out->solver_package=PASO_PASO;
50	out->solver=NULL;
51	out->val=NULL;
52	out->trilinos_data=NULL;
53	out->reference_counter=1;
54	out->type=type;
55	/* ====== compressed sparse columns === */
56	if (type & MATRIX_FORMAT_CSC) {
57	if (type & MATRIX_FORMAT_SYM) {
58	Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
59	return NULL;
60	} else {
61	if ((type & MATRIX_FORMAT_BLK1) \|\| row_block_size!=col_block_size \|\| col_block_size>3) {
62	if (type & MATRIX_FORMAT_OFFSET1) {
63	out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,col_block_size,row_block_size);
64	} else {
65	out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,col_block_size,row_block_size);
66	}
67	out->row_block_size=1;
68	out->col_block_size=1;
69	} else {
70	if ( (type & MATRIX_FORMAT_OFFSET1) ==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
71	out->pattern=Paso_SystemMatrixPattern_reference(pattern);
72	} else {
73	out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT,1,1);
74	}
75	out->row_block_size=row_block_size;
76	out->col_block_size=col_block_size;
77	}
78	}
79	out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
80	out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
81	} else {
82	/* ====== compressed sparse row === */
83	if (type & MATRIX_FORMAT_SYM) {
84	Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
85	return NULL;
86	} else {
87	if ((type & MATRIX_FORMAT_BLK1) \|\| row_block_size!=col_block_size \|\| col_block_size>3) {
88	if (type & MATRIX_FORMAT_OFFSET1) {
89	out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,row_block_size,col_block_size);
90	} else {
91	out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,row_block_size,col_block_size);
92	}
93	out->row_block_size=1;
94	out->col_block_size=1;
95	} else {
96	if ((type & MATRIX_FORMAT_OFFSET1)==(pattern->type & PATTERN_FORMAT_OFFSET1)) {
97	out->pattern=Paso_SystemMatrixPattern_reference(pattern);
98	} else {
99	out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT,1,1);
100	}
101	out->row_block_size=row_block_size;
102	out->col_block_size=col_block_size;
103	}
104	}
105	out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
106	out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
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->myLen=(size_t)(out->pattern->myLen)*(size_t)(out->block_size);
113
114	out->numRows = out->row_distribution->myNumComponents;
115	out->myNumRows = out->row_distribution->numComponents;
116	out->numCols = out->col_distribution->myNumComponents;
117	out->myNumCols = out->col_distribution->numComponents;
118	out->mpi_info = Paso_MPIInfo_getReference(out->pattern->mpi_info);
119	/* allocate memory for matrix entries */
120	if (type & MATRIX_FORMAT_TRILINOS_CRS) {
121	Paso_TRILINOS_alloc(out->trilinos_data, out->pattern,out->row_block_size,out->col_block_size);
122	} else {
123	if (type & MATRIX_FORMAT_CSC) {
124	n_norm = out->myNumCols * out->col_block_size;
125	} else {
126	n_norm = out->myNumRows * out->row_block_size;
127	}
128	out->val=MEMALLOC(out->myLen,double);
129	out->normalizer=MEMALLOC(n_norm,double);
130	out->normalizer_is_valid=FALSE;
131	if (! Paso_checkPtr(out->val)) Paso_SystemMatrix_setValues(out,DBLE(0));
132	if (! Paso_checkPtr(out->normalizer)) {
133	#pragma omp parallel for private(i) schedule(static)
134	for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
135	}
136	}
137	}
138	/* all done: */
139	if (! Paso_noError()) {
140	Paso_SystemMatrix_dealloc(out);
141	return NULL;
142	} else {
143	#ifdef Paso_TRACE
144	printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
145	printf("Paso_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->numRows,(long)out->numCols);
146	#endif
147	return out;
148	}
149	}
150
151	/* returns a reference to Paso_SystemMatrix in */
152
153	Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
154	if (in!=NULL) ++(in->reference_counter);
155	return NULL;
156	}
157
158	/* deallocates a SystemMatrix: */
159
160	void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {
161	if (in!=NULL) {
162	in->reference_counter--;
163	if (in->reference_counter<=0) {
164	Paso_SystemMatrixPattern_dealloc(in->pattern);
165	Paso_Distribution_free(in->row_distribution);
166	Paso_Distribution_free(in->col_distribution);
167	Paso_MPIInfo_dealloc(in->mpi_info);
168	MEMFREE(in->val);
169	MEMFREE(in->normalizer);
170	Paso_TRILINOS_free(in->trilinos_data);
171	Paso_solve_free(in);
172	MEMFREE(in);
173	#ifdef Paso_TRACE
174	printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");
175	#endif
176	}
177	}
178	}
179	/*
180	* $Log$
181	* Revision 1.2 2005/09/15 03:44:38 jgs
182	* Merge of development branch dev-02 back to main trunk on 2005-09-15
183	*
184	* Revision 1.1.2.2 2005/09/07 00:59:08 gross
185	* some inconsistent renaming fixed to make the linking work.
186	*
187	* Revision 1.1.2.1 2005/09/05 06:29:47 gross
188	* These files have been extracted from finley to define a stand alone libray for iterative
189	* linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
190	* has not been tested yet.
191	*
192	*
193	*/
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision