sundials_direct.h 13 KB

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