00001 /* lah_block.h - local header for matrix storage and operations package */ 00002 /* based on block storage and providing basic operations for */ 00003 /* standard iterative methods (Jacobi, Gauss-Seidel, additive Schwarz, */ 00004 /* multiplicative Schwarz) */ 00005 /* Assumption: Each matrix is stored with the corresponding preconditioner */ 00006 00007 00008 #ifndef _lah_block_ 00009 #define _lah_block_ 00010 00011 /*** External Parameters ***/ 00012 00013 #define LAC_MAX_MATRICES 20 00014 00015 00016 #ifdef __cplusplus 00017 extern "C" { 00018 #endif 00019 00020 /*** Data types ***/ 00021 00022 /* definition of type itt_block */ 00023 typedef struct { 00024 /* number of degrees of freedom in a single diagonal elementary block */ 00025 int Ndof; 00026 /* list of neighbor blocks, Lngb[0] - number of neighbors */ 00027 int *Lngb; 00028 /* position of the first dof in the global rhs vector */ 00029 int Posg; 00030 /* array for pivoting information*/ 00031 int *Ips; 00032 /* stiffness matrix for diagonal elementary block */ 00033 double *Dia; 00034 /* RHS vector for diagonal elementary block*/ 00035 double *Rhs; 00036 /* stiffness matrices for off diagonal elementary blocks */ 00037 double **Aux; 00038 } itt_blocks; 00039 00040 /* definition of type itt_blocks_dia */ 00041 typedef struct { 00042 /* list of elementary blocks' numbers, Lsmall[0] - number of blocks */ 00043 int *Lsmall; 00044 /* list of elements in the subdoamin, Lelem[0] - number of elements */ 00045 int *Lelem; 00046 /* list of neighbors of the subdoamin, Lneig[0] - number of neighbors */ 00047 int *Lneig; 00048 /* list of positions of first dofs from elementary blocks in diagonal 00049 subarray; Lpos[0] - total number of dofs for dia */ 00050 int *Lpos; 00051 /* list of lower subdiagonal neighbors - for ILU preconditioning, */ 00052 /* Llowerneig[0] - number of neighbors */ 00053 int *Llowerneig; 00054 /* list of upper subdiagonal neighbors - for ILU preconditioning, */ 00055 /* Lupperneig[0] - number of neighbors */ 00056 int *Lupperneig; 00057 /* array for pivoting information*/ 00058 int *Ips; 00059 /* factorized diagonal blocks (subarray) */ 00060 double *Dia; 00061 /* factorized off diagonal blocks */ 00062 double **Aux; 00063 } itt_blocks_dia; 00064 00065 00066 typedef struct { 00067 00068 /* control variables */ 00069 int Block_type; /* block types: number of nodes in a block */ 00070 /* or some other indicator application dependent */ 00071 int Precon; /* type of preconditioner determining the storage format */ 00072 /* of the preconditioner matrix Block_dia */ 00073 int ILU_k; // number of neighborhood rings for ILU(k) 00074 /* parameters */ 00075 int Max_SM_size; /* the maximal number of dofs in a stiffness matrix */ 00076 int Nrblocks; /* total number of small blocks */ 00077 int Nrblocks_dia; /* total number of diagonal blocks */ 00078 int Nrdofgl; /* total number of degrees of freedom */ 00079 00080 /* blocks for storing system array and preconditioner arrays */ 00081 itt_blocks **Block; /* array of pointers to small blocks */ 00082 itt_blocks_dia **Block_dia; /* array of pointers to diagonal blocks */ 00083 00084 } itt_matrices; 00085 00086 00087 /* GLOBAL VARIABLES */ 00088 extern int itv_nr_matrices; /* the number of solvers in the problem */ 00089 extern int itv_cur_matrix_id; /* ID of the current problem */ 00090 extern itt_matrices itv_matrices[LAC_MAX_MATRICES]; /* array of solvers */ 00091 00092 00093 /* INTERFACE IMPLEMENTATIONS - the same routines as in lah_intf.h but for _block */ 00094 00098 int lar_allocate_SM_and_LV_block( // returns: matrix index in itv_matrices array 00099 int Nrdof_glob, /* in: total number of DOFs */ 00100 int Max_SM_size, /* maximal size of element stiffness matrix */ 00101 int Nrblocks, /* in: number of DOF blocks */ 00102 int* Nrdofbl, /* in: list of numbers of dofs in a block */ 00103 int* Posglob, /* in: list of global numbers of first dof */ 00104 int* Nroffbl, /* in: list of numbers of off diagonal blocks */ 00105 // if Nroffbl[iblock]<=0 - iblock is a ghost block with no rows associated 00106 int** L_offbl /* in: list of lists of off diagonal blocks */ 00107 ); 00108 00112 extern int lar_initialize_SM_and_LV_block( 00113 int Matrix_id, /* in: matrix ID */ 00114 int Scope /* in: indicator for the scope of computations: */ 00115 /* LAC_SCOPE_SM_AND_LV */ 00116 /* LAC_SCOPE_LV - do not touch SM! */ 00117 ); 00118 00122 extern double lar_get_storage_block( /* returns: storage in MB */ 00123 int Matrix_id /* in: matrix ID */ 00124 ); 00125 00131 extern int lar_assemble_SM_and_LV_block( 00132 /* returns: >=0 - success code, <0 - error code */ 00133 int Matrix_id, /* in: matrix ID */ 00134 int Scope, /* in: indicator for the scope of computations: */ 00135 /* LAC_SCOPE_SM_AND_LV */ 00136 /* LAC_SCOPE_LV - do not touch SM! */ 00137 int Nr_dof_bl, /* in: number of global dof blocks */ 00138 /* associated with the local stiffness matrix */ 00139 int* L_bl_id, /* in: list of dof blocks' IDs */ 00140 int* L_bl_nrdof, /* in: list of blocks' numbers of dof */ 00141 double* Stiff_mat, /* in: stiffness matrix stored columnwise */ 00142 double* Rhs_vect, /* in: rhs vector */ 00143 char* Rewr_dofs /* in: flag to rewrite or sum up entries */ 00144 /* 'T' - true, rewrite entries when assembling */ 00145 /* 'F' - false, sum up entries when assembling */ 00146 ); 00147 00148 00152 extern int lar_allocate_preconditioner_block( 00153 /* returns: >0 number of diagonal blocks */ 00154 /* <=0 - error */ 00155 int Matrix_id, /* in: matrix ID */ 00156 int Precon, /* in: type of preconditioner (lah_block.h line circa 45) */ 00157 int Block_type /* in: number of elementary DOF blocks in a solver block */ 00158 ); 00159 00163 extern int lar_fill_preconditioner_block( 00164 int Matrix_id /* in: matrix ID */ 00165 ); 00166 00170 extern int lar_free_preconditioner_block( 00171 int Matrix_id /* in: matrix ID */ 00172 ); 00173 00177 extern int lar_free_SM_and_LV_block( 00178 int Matrix_id /* in: matrix ID */ 00179 ); 00180 00181 00186 extern void lar_compute_residual_block ( 00187 int Matrix_id, /* in: matrix ID */ 00188 int Use_rhs, /* in: indicator whether to use RHS */ 00189 int Ini_zero, /* in: flag for zero initial guess */ 00190 int Ndof, /* in: number of unknowns (components of x) */ 00191 double* X, /* in: initial guess vector */ 00192 double* B, /* in: the rhs vector, if NULL take rhs */ 00193 /* from block data structure */ 00194 double* V /* out: residual, v = ( b-Ax ) */ 00195 ); 00196 00202 extern void lar_compute_preconditioned_residual_block ( 00203 int Matrix_id, /* in: matrix ID */ 00204 int Use_rhs, /* in: indicator whether to use RHS */ 00205 int Ini_zero, /* in: flag for zero initial guess */ 00206 int Ndof, /* in: number of unknowns (components of x) */ 00207 double* X, /* in: initial guess vector */ 00208 double* B, /* in: the rhs vector, if NULL take rhs */ 00209 /* from block data structure */ 00210 double* V /* out: preconditioned residual, v = M^-1*(b-Ax) */ 00211 ); 00212 00218 extern void lar_perform_BJ_or_GS_iterations_block( 00219 int Matrix_id, /* in: matrix ID */ 00220 int Use_rhs, /* in: 0 - no rhs, 1 - with rhs */ 00221 int Ini_zero, /* in: flag for zero initial guess */ 00222 int Nr_prec, /* in: number of preconditioner iterations */ 00223 int Ndof, /* in: number of unknowns (components of v*) */ 00224 double* V, /* in,out: vector of unknowns updated */ 00225 /* during the loop over subdomains */ 00226 double* B /* in: the rhs vector, if NULL take rhs */ 00227 /* from block data structure */ 00228 ); 00229 00230 00235 extern void lar_perform_rhsub_block( 00236 int Matrix_id, /* in: matrix ID */ 00237 int Ndof, /* in: number of unknowns (components of v*) */ 00238 double* V, /* in,out: vector of unknowns updated */ 00239 /* during the loop over subdomains */ 00240 double* B /* in: the rhs vector, if NULL take rhs */ 00241 /* from block data structure */ 00242 ); 00243 00244 00245 00246 // internal utilities 00247 /* 00248 lar_util_dvector - to allocate a double vector: name[0..ncom-1]: 00249 lar_util_ivector - to allocate an integer vector: name[0..ncom-1]: 00250 lar_util_imatrix - to allocate an integer matrix name[0..nrow-1][0..ncol-1]: 00251 name=imatrix(nrow,ncol,error_text) 00252 lar_util_dmatrix - to allocate a double matrix name[0..nrow-1][0..ncol-1]: 00253 name=imatrix(nrow,ncol,error_text) 00254 lar_util_chk_list - to check whether a number is on the list 00255 lar_util_put_list - to put Num on the list List with length Ll 00256 lar_util_d_zero - to zero a double vector 00257 lar_util_i_zero - to zero an integer vector 00258 lar_util_sort - to heap-sort an array 00259 lar_util_dgetrf - quasi-LU decomposition of a matrix 00260 lar_util_dgetrs - to perform forward reduction and back substitution 00261 of the RHS vector for solving a system of linear equations 00262 */ 00263 00264 00265 extern double *lar_util_dvector( 00266 /* return: pointer to allocated vector */ 00267 int Ncom, /* in: number of components */ 00268 char Error_text[]/* in: error text to be printed */ 00269 ); 00270 00271 extern int *lar_util_ivector( 00272 /* return: pointer to allocated vector */ 00273 int Ncom, /* in: number of components */ 00274 char Error_text[]/* in: error text to be printed */ 00275 ); 00276 00277 extern int **lar_util_imatrix( /* returns: pointer to array of pointers to integers */ 00278 int Nrow, /* in: number of rows */ 00279 int Ncol, /* in: number of columns */ 00280 char Error_text[]/* in: text to print in case of error */ 00281 ); 00282 00283 extern double **lar_util_dmatrix( /* returns: pointer to array of pointers to doubles */ 00284 int Nrow, /* in: number of rows */ 00285 int Ncol, /* in: number of columns */ 00286 char Error_text[]/* in: text to print in case of error */ 00287 ); 00288 00289 extern int lar_util_chk_list( /* returns: */ 00290 /* >0 - position on the list */ 00291 /* 0 - not found on the list */ 00292 int Num, /* number to be checked */ 00293 int* List, /* list of numbers */ 00294 int Ll /* length of the list */ 00295 ); 00296 00297 extern int lar_util_put_list( /* returns*/ 00298 /* >0 - position already occupied on the list */ 00299 /* 0 - put on the list */ 00300 /* -1 - list full, not found on the list */ 00301 int Num, /* in: number to put on the list */ 00302 int* List, /* in: list */ 00303 int Ll /* in: total list's lengths */ 00304 ); 00305 00306 extern void lar_util_d_zero(double *Vec, int Num); 00307 00308 extern void lar_util_i_zero(int *Vec, int Num); 00309 00310 extern void lar_util_sort( 00311 int *Ind_array, /* in/out: index array for sorting */ 00312 double *Val_array /* in: array of values used for sorting */ 00313 ); 00314 00315 extern void lar_util_dgetrf(double* a, int m, int* ips); 00316 00317 extern void lar_util_dgetrs(double* a, int m, double* b, double* x, int* ips); 00318 00319 00320 #ifdef __cplusplus 00321 } 00322 #endif 00323 00324 #endif