[16222] | 1 | /* ----------------------------------------------------------------- |
---|
| 2 | * Programmer: Radu Serban @ LLNL |
---|
| 3 | * ----------------------------------------------------------------- |
---|
| 4 | * LLNS Copyright Start |
---|
| 5 | * Copyright (c) 2014, Lawrence Livermore National Security |
---|
| 6 | * This work was performed under the auspices of the U.S. Department |
---|
| 7 | * of Energy by Lawrence Livermore National Laboratory in part under |
---|
| 8 | * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344. |
---|
| 9 | * Produced at the Lawrence Livermore National Laboratory. |
---|
| 10 | * All rights reserved. |
---|
| 11 | * For details, see the LICENSE file. |
---|
| 12 | * LLNS 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, sunindextype ml, sunindextype smu); |
---|
| 221 | |
---|
| 222 | /* |
---|
| 223 | * ----------------------------------------------------------------- |
---|
| 224 | * Functions: DestroyMat |
---|
| 225 | * ----------------------------------------------------------------- |
---|
| 226 | * DestroyMat frees the memory allocated by NewDenseMat or NewBandMat |
---|
| 227 | * ----------------------------------------------------------------- |
---|
| 228 | */ |
---|
| 229 | |
---|
| 230 | SUNDIALS_EXPORT void DestroyMat(DlsMat A); |
---|
| 231 | |
---|
| 232 | /* |
---|
| 233 | * ----------------------------------------------------------------- |
---|
| 234 | * Function: NewIntArray |
---|
| 235 | * ----------------------------------------------------------------- |
---|
| 236 | * NewIntArray allocates memory an array of N int's and returns |
---|
| 237 | * the pointer to the memory it allocates. If the request for |
---|
| 238 | * memory storage cannot be satisfied, it returns NULL. |
---|
| 239 | * ----------------------------------------------------------------- |
---|
| 240 | */ |
---|
| 241 | |
---|
| 242 | SUNDIALS_EXPORT int *NewIntArray(int N); |
---|
| 243 | |
---|
| 244 | /* |
---|
| 245 | * ----------------------------------------------------------------- |
---|
| 246 | * Function: NewIndexArray |
---|
| 247 | * ----------------------------------------------------------------- |
---|
| 248 | * NewIndexArray allocates memory an array of N sunindextype's and |
---|
| 249 | * returns the pointer to the memory it allocates. If the request |
---|
| 250 | * for memory storage cannot be satisfied, it returns NULL. |
---|
| 251 | * ----------------------------------------------------------------- |
---|
| 252 | */ |
---|
| 253 | |
---|
| 254 | SUNDIALS_EXPORT sunindextype *NewIndexArray(sunindextype N); |
---|
| 255 | |
---|
| 256 | /* |
---|
| 257 | * ----------------------------------------------------------------- |
---|
| 258 | * Function: NewRealArray |
---|
| 259 | * ----------------------------------------------------------------- |
---|
| 260 | * NewRealArray allocates memory an array of N realtype and returns |
---|
| 261 | * the pointer to the memory it allocates. If the request for |
---|
| 262 | * memory storage cannot be satisfied, it returns NULL. |
---|
| 263 | * ----------------------------------------------------------------- |
---|
| 264 | */ |
---|
| 265 | |
---|
| 266 | SUNDIALS_EXPORT realtype *NewRealArray(sunindextype N); |
---|
| 267 | |
---|
| 268 | /* |
---|
| 269 | * ----------------------------------------------------------------- |
---|
| 270 | * Function: DestroyArray |
---|
| 271 | * ----------------------------------------------------------------- |
---|
| 272 | * DestroyArray frees memory allocated by NewIntArray, NewIndexArray, |
---|
| 273 | * or NewRealArray. |
---|
| 274 | * ----------------------------------------------------------------- |
---|
| 275 | */ |
---|
| 276 | |
---|
| 277 | SUNDIALS_EXPORT void DestroyArray(void *p); |
---|
| 278 | |
---|
| 279 | /* |
---|
| 280 | * ----------------------------------------------------------------- |
---|
| 281 | * Function : AddIdentity |
---|
| 282 | * ----------------------------------------------------------------- |
---|
| 283 | * AddIdentity adds 1.0 to the main diagonal (A_ii, i=0,1,...,N-1) of |
---|
| 284 | * the M-by-N matrix A (M>= N) and stores the result back in A. |
---|
| 285 | * AddIdentity is typically used with square matrices. |
---|
| 286 | * AddIdentity does not check for M >= N and therefore a segmentation |
---|
| 287 | * fault will occur if M < N! |
---|
| 288 | * ----------------------------------------------------------------- |
---|
| 289 | */ |
---|
| 290 | |
---|
| 291 | SUNDIALS_EXPORT void AddIdentity(DlsMat A); |
---|
| 292 | |
---|
| 293 | /* |
---|
| 294 | * ----------------------------------------------------------------- |
---|
| 295 | * Function : SetToZero |
---|
| 296 | * ----------------------------------------------------------------- |
---|
| 297 | * SetToZero sets all the elements of the M-by-N matrix A to 0.0. |
---|
| 298 | * ----------------------------------------------------------------- |
---|
| 299 | */ |
---|
| 300 | |
---|
| 301 | SUNDIALS_EXPORT void SetToZero(DlsMat A); |
---|
| 302 | |
---|
| 303 | /* |
---|
| 304 | * ----------------------------------------------------------------- |
---|
| 305 | * Functions: PrintMat |
---|
| 306 | * ----------------------------------------------------------------- |
---|
| 307 | * This function prints the M-by-N (dense or band) matrix A to |
---|
| 308 | * outfile as it would normally appear on paper. |
---|
| 309 | * It is intended as debugging tools with small values of M and N. |
---|
| 310 | * The elements are printed using the %g/%lg/%Lg option. |
---|
| 311 | * A blank line is printed before and after the matrix. |
---|
| 312 | * ----------------------------------------------------------------- |
---|
| 313 | */ |
---|
| 314 | |
---|
| 315 | SUNDIALS_EXPORT void PrintMat(DlsMat A, FILE *outfile); |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | /* |
---|
| 319 | * ================================================================== |
---|
| 320 | * Exported function prototypes (functions working on realtype**) |
---|
| 321 | * ================================================================== |
---|
| 322 | */ |
---|
| 323 | |
---|
| 324 | SUNDIALS_EXPORT realtype **newDenseMat(sunindextype m, sunindextype n); |
---|
| 325 | SUNDIALS_EXPORT realtype **newBandMat(sunindextype n, sunindextype smu, sunindextype ml); |
---|
| 326 | SUNDIALS_EXPORT void destroyMat(realtype **a); |
---|
| 327 | SUNDIALS_EXPORT int *newIntArray(int n); |
---|
| 328 | SUNDIALS_EXPORT sunindextype *newIndexArray(sunindextype n); |
---|
| 329 | SUNDIALS_EXPORT realtype *newRealArray(sunindextype m); |
---|
| 330 | SUNDIALS_EXPORT void destroyArray(void *v); |
---|
| 331 | |
---|
| 332 | |
---|
| 333 | #ifdef __cplusplus |
---|
| 334 | } |
---|
| 335 | #endif |
---|
| 336 | |
---|
| 337 | #endif |
---|