1 /* ----------------------------------------------------------------- 2 * Programmer: Radu Serban @ LLNL 3 * ----------------------------------------------------------------- 4 * SUNDIALS Copyright Start 5 * Copyright (c) 2002-2021, Lawrence Livermore National Security 6 * and Southern Methodist University. 7 * All rights reserved. 8 * 9 * See the top-level LICENSE and NOTICE files for details. 10 * 11 * SPDX-License-Identifier: BSD-3-Clause 12 * SUNDIALS Copyright End 13 * ----------------------------------------------------------------- 14 * This header file contains definitions and declarations for use by 15 * generic direct linear solvers for Ax = b. It defines types for 16 * dense and banded matrices and corresponding accessor macros. 17 * -----------------------------------------------------------------*/ 18 19 #ifndef _SUNDIALS_DIRECT_H 20 #define _SUNDIALS_DIRECT_H 21 22 #include <stdio.h> 23 #include <sundials/sundials_types.h> 24 25 #ifdef __cplusplus /* wrapper to enable C++ usage */ 26 extern "C" { 27 #endif 28 29 /* 30 * ================================================================= 31 * C O N S T A N T S 32 * ================================================================= 33 */ 34 35 /* 36 * SUNDIALS_DENSE: dense matrix 37 * SUNDIALS_BAND: banded matrix 38 */ 39 40 #define SUNDIALS_DENSE 1 41 #define SUNDIALS_BAND 2 42 43 /* 44 * ================================================================== 45 * Type definitions 46 * ================================================================== 47 */ 48 49 /* 50 * ----------------------------------------------------------------- 51 * Type : DlsMat 52 * ----------------------------------------------------------------- 53 * The type DlsMat is defined to be a pointer to a structure 54 * with various sizes, a data field, and an array of pointers to 55 * the columns which defines a dense or band matrix for use in 56 * direct linear solvers. The M and N fields indicates the number 57 * of rows and columns, respectively. The data field is a one 58 * dimensional array used for component storage. The cols field 59 * stores the pointers in data for the beginning of each column. 60 * ----------------------------------------------------------------- 61 * For DENSE matrices, the relevant fields in DlsMat are: 62 * type = SUNDIALS_DENSE 63 * M - number of rows 64 * N - number of columns 65 * ldim - leading dimension (ldim >= M) 66 * data - pointer to a contiguous block of realtype variables 67 * ldata - length of the data array =ldim*N 68 * cols - array of pointers. cols[j] points to the first element 69 * of the j-th column of the matrix in the array data. 70 * 71 * The elements of a dense matrix are stored columnwise (i.e. columns 72 * are stored one on top of the other in memory). 73 * If A is of type DlsMat, then the (i,j)th element of A (with 74 * 0 <= i < M and 0 <= j < N) is given by (A->data)[j*n+i]. 75 * 76 * The DENSE_COL and DENSE_ELEM macros below allow a user to access 77 * efficiently individual matrix elements without writing out explicit 78 * data structure references and without knowing too much about the 79 * underlying element storage. The only storage assumption needed is 80 * that elements are stored columnwise and that a pointer to the 81 * jth column of elements can be obtained via the DENSE_COL macro. 82 * ----------------------------------------------------------------- 83 * For BAND matrices, the relevant fields in DlsMat are: 84 * type = SUNDIALS_BAND 85 * M - number of rows 86 * N - number of columns 87 * mu - upper bandwidth, 0 <= mu <= min(M,N) 88 * ml - lower bandwidth, 0 <= ml <= min(M,N) 89 * s_mu - storage upper bandwidth, mu <= s_mu <= N-1. 90 * The dgbtrf routine writes the LU factors into the storage 91 * for A. The upper triangular factor U, however, may have 92 * an upper bandwidth as big as MIN(N-1,mu+ml) because of 93 * partial pivoting. The s_mu field holds the upper 94 * bandwidth allocated for A. 95 * ldim - leading dimension (ldim >= s_mu) 96 * data - pointer to a contiguous block of realtype variables 97 * ldata - length of the data array =ldim*(s_mu+ml+1) 98 * cols - array of pointers. cols[j] points to the first element 99 * of the j-th column of the matrix in the array data. 100 * 101 * The BAND_COL, BAND_COL_ELEM, and BAND_ELEM macros below allow a 102 * user to access individual matrix elements without writing out 103 * explicit data structure references and without knowing too much 104 * about the underlying element storage. The only storage assumption 105 * needed is that elements are stored columnwise and that a pointer 106 * into the jth column of elements can be obtained via the BAND_COL 107 * macro. The BAND_COL_ELEM macro selects an element from a column 108 * which has already been isolated via BAND_COL. The macro 109 * BAND_COL_ELEM allows the user to avoid the translation 110 * from the matrix location (i,j) to the index in the array returned 111 * by BAND_COL at which the (i,j)th element is stored. 112 * ----------------------------------------------------------------- 113 */ 114 115 typedef struct _DlsMat { 116 int type; 117 sunindextype M; 118 sunindextype N; 119 sunindextype ldim; 120 sunindextype mu; 121 sunindextype ml; 122 sunindextype s_mu; 123 realtype *data; 124 sunindextype ldata; 125 realtype **cols; 126 } *DlsMat; 127 128 /* 129 * ================================================================== 130 * Data accessor macros 131 * ================================================================== 132 */ 133 134 /* 135 * ----------------------------------------------------------------- 136 * DENSE_COL and DENSE_ELEM 137 * ----------------------------------------------------------------- 138 * 139 * DENSE_COL(A,j) references the jth column of the M-by-N dense 140 * matrix A, 0 <= j < N. The type of the expression DENSE_COL(A,j) 141 * is (realtype *). After the assignment col_j = DENSE_COL(A,j), 142 * col_j may be treated as an array indexed from 0 to M-1. 143 * The (i,j)-th element of A is thus referenced by col_j[i]. 144 * 145 * DENSE_ELEM(A,i,j) references the (i,j)th element of the dense 146 * M-by-N matrix A, 0 <= i < M ; 0 <= j < N. 147 * 148 * ----------------------------------------------------------------- 149 */ 150 151 #define DENSE_COL(A,j) ((A->cols)[j]) 152 #define DENSE_ELEM(A,i,j) ((A->cols)[j][i]) 153 154 /* 155 * ----------------------------------------------------------------- 156 * BAND_COL, BAND_COL_ELEM, and BAND_ELEM 157 * ----------------------------------------------------------------- 158 * 159 * BAND_COL(A,j) references the diagonal element of the jth column 160 * of the N by N band matrix A, 0 <= j <= N-1. The type of the 161 * expression BAND_COL(A,j) is realtype *. The pointer returned by 162 * the call BAND_COL(A,j) can be treated as an array which is 163 * indexed from -(A->mu) to (A->ml). 164 * 165 * BAND_COL_ELEM references the (i,j)th entry of the band matrix A 166 * when used in conjunction with BAND_COL. The index (i,j) should 167 * satisfy j-(A->mu) <= i <= j+(A->ml). 168 * 169 * BAND_ELEM(A,i,j) references the (i,j)th element of the M-by-N 170 * band matrix A, where 0 <= i,j <= N-1. The location (i,j) should 171 * further satisfy j-(A->mu) <= i <= j+(A->ml). 172 * 173 * ----------------------------------------------------------------- 174 */ 175 176 #define BAND_COL(A,j) (((A->cols)[j])+(A->s_mu)) 177 #define BAND_COL_ELEM(col_j,i,j) (col_j[(i)-(j)]) 178 #define BAND_ELEM(A,i,j) ((A->cols)[j][(i)-(j)+(A->s_mu)]) 179 180 /* 181 * ================================================================== 182 * Exported function prototypes (functions working on dlsMat) 183 * ================================================================== 184 */ 185 186 /* 187 * ----------------------------------------------------------------- 188 * Function: NewDenseMat 189 * ----------------------------------------------------------------- 190 * NewDenseMat allocates memory for an M-by-N dense matrix and 191 * returns the storage allocated (type DlsMat). NewDenseMat 192 * returns NULL if the request for matrix storage cannot be 193 * satisfied. See the above documentation for the type DlsMat 194 * for matrix storage details. 195 * ----------------------------------------------------------------- 196 */ 197 198 SUNDIALS_EXPORT DlsMat NewDenseMat(sunindextype M, sunindextype N); 199 200 /* 201 * ----------------------------------------------------------------- 202 * Function: NewBandMat 203 * ----------------------------------------------------------------- 204 * NewBandMat allocates memory for an M-by-N band matrix 205 * with upper bandwidth mu, lower bandwidth ml, and storage upper 206 * bandwidth smu. Pass smu as follows depending on whether A will 207 * be LU factored: 208 * 209 * (1) Pass smu = mu if A will not be factored. 210 * 211 * (2) Pass smu = MIN(N-1,mu+ml) if A will be factored. 212 * 213 * NewBandMat returns the storage allocated (type DlsMat) or 214 * NULL if the request for matrix storage cannot be satisfied. 215 * See the documentation for the type DlsMat for matrix storage 216 * details. 217 * ----------------------------------------------------------------- 218 */ 219 220 SUNDIALS_EXPORT DlsMat NewBandMat(sunindextype N, sunindextype mu, 221 sunindextype ml, sunindextype smu); 222 223 /* 224 * ----------------------------------------------------------------- 225 * Functions: DestroyMat 226 * ----------------------------------------------------------------- 227 * DestroyMat frees the memory allocated by NewDenseMat or NewBandMat 228 * ----------------------------------------------------------------- 229 */ 230 231 SUNDIALS_EXPORT void DestroyMat(DlsMat A); 232 233 /* 234 * ----------------------------------------------------------------- 235 * Function: NewIntArray 236 * ----------------------------------------------------------------- 237 * NewIntArray allocates memory an array of N int's and returns 238 * the pointer to the memory it allocates. If the request for 239 * memory storage cannot be satisfied, it returns NULL. 240 * ----------------------------------------------------------------- 241 */ 242 243 SUNDIALS_EXPORT int *NewIntArray(int N); 244 245 /* 246 * ----------------------------------------------------------------- 247 * Function: NewIndexArray 248 * ----------------------------------------------------------------- 249 * NewIndexArray allocates memory an array of N sunindextype's and 250 * returns the pointer to the memory it allocates. If the request 251 * for memory storage cannot be satisfied, it returns NULL. 252 * ----------------------------------------------------------------- 253 */ 254 255 SUNDIALS_EXPORT sunindextype *NewIndexArray(sunindextype N); 256 257 /* 258 * ----------------------------------------------------------------- 259 * Function: NewRealArray 260 * ----------------------------------------------------------------- 261 * NewRealArray allocates memory an array of N realtype and returns 262 * the pointer to the memory it allocates. If the request for 263 * memory storage cannot be satisfied, it returns NULL. 264 * ----------------------------------------------------------------- 265 */ 266 267 SUNDIALS_EXPORT realtype *NewRealArray(sunindextype N); 268 269 /* 270 * ----------------------------------------------------------------- 271 * Function: DestroyArray 272 * ----------------------------------------------------------------- 273 * DestroyArray frees memory allocated by NewIntArray, NewIndexArray, 274 * or NewRealArray. 275 * ----------------------------------------------------------------- 276 */ 277 278 SUNDIALS_EXPORT void DestroyArray(void *p); 279 280 /* 281 * ----------------------------------------------------------------- 282 * Function : AddIdentity 283 * ----------------------------------------------------------------- 284 * AddIdentity adds 1.0 to the main diagonal (A_ii, i=0,1,...,N-1) of 285 * the M-by-N matrix A (M>= N) and stores the result back in A. 286 * AddIdentity is typically used with square matrices. 287 * AddIdentity does not check for M >= N and therefore a segmentation 288 * fault will occur if M < N! 289 * ----------------------------------------------------------------- 290 */ 291 292 SUNDIALS_EXPORT void AddIdentity(DlsMat A); 293 294 /* 295 * ----------------------------------------------------------------- 296 * Function : SetToZero 297 * ----------------------------------------------------------------- 298 * SetToZero sets all the elements of the M-by-N matrix A to 0.0. 299 * ----------------------------------------------------------------- 300 */ 301 302 SUNDIALS_EXPORT void SetToZero(DlsMat A); 303 304 /* 305 * ----------------------------------------------------------------- 306 * Functions: PrintMat 307 * ----------------------------------------------------------------- 308 * This function prints the M-by-N (dense or band) matrix A to 309 * outfile as it would normally appear on paper. 310 * It is intended as debugging tools with small values of M and N. 311 * The elements are printed using the %g/%lg/%Lg option. 312 * A blank line is printed before and after the matrix. 313 * ----------------------------------------------------------------- 314 */ 315 316 SUNDIALS_EXPORT void PrintMat(DlsMat A, FILE *outfile); 317 318 319 /* 320 * ================================================================== 321 * Exported function prototypes (functions working on realtype**) 322 * ================================================================== 323 */ 324 325 SUNDIALS_EXPORT realtype **newDenseMat(sunindextype m, sunindextype n); 326 SUNDIALS_EXPORT realtype **newBandMat(sunindextype n, sunindextype smu, 327 sunindextype ml); 328 SUNDIALS_EXPORT void destroyMat(realtype **a); 329 SUNDIALS_EXPORT int *newIntArray(int n); 330 SUNDIALS_EXPORT sunindextype *newIndexArray(sunindextype n); 331 SUNDIALS_EXPORT realtype *newRealArray(sunindextype m); 332 SUNDIALS_EXPORT void destroyArray(void *v); 333 334 335 #ifdef __cplusplus 336 } 337 #endif 338 339 #endif 340