Ice Sheet System Model  4.18
Code documentation
Functions
gslincludes.h File Reference
#include "../../shared/Numerics/types.h"

Go to the source code of this file.

Functions

void DenseGslSolve (IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B, int n)
 
void DenseGslSolve (IssmDouble **pX, IssmDouble *Kff, int Kff_M, int Kff_N, IssmDouble *pf, int pf_M, Parameters *parameters)
 
void SolverxSeq (IssmPDouble *X, IssmPDouble *A, IssmPDouble *B, int n)
 

Function Documentation

◆ DenseGslSolve() [1/2]

void DenseGslSolve ( IssmPDouble **  pX,
IssmPDouble A,
IssmPDouble B,
int  n 
)

Definition at line 25 of file DenseGslSolve.cpp.

25  { /*{{{*/
26 
27  /*Intermediary: */
28  IssmPDouble *X = xNew<IssmPDouble>(n);
29  SolverxSeq(X,A,B,n);
30 
31  /*allocate output pointers: */
32  *pX=X;
33 }

◆ DenseGslSolve() [2/2]

void DenseGslSolve ( IssmDouble **  pX,
IssmDouble Kff,
int  Kff_M,
int  Kff_N,
IssmDouble pf,
int  pf_M,
Parameters parameters 
)

Definition at line 35 of file DenseGslSolve.cpp.

35  { /*{{{*/
36 
37  /*Intermediary: */
38 
39  if(Kff_N!=pf_M)_error_("Right hand side vector of size " << pf_M << ", when matrix is of size " << Kff_M << "-" << Kff_N << " !");
40  if(Kff_M!=Kff_N)_error_("Stiffness matrix should be square!");
41 
42  IssmPDouble *x = xNew<IssmPDouble>(Kff_N);
43  SolverxSeq(x,Kff,pf,Kff_N);
44 
45  /*allocate output pointers: */
46  *px=x;
47 }

◆ SolverxSeq()

void SolverxSeq ( IssmPDouble X,
IssmPDouble A,
IssmPDouble B,
int  n 
)

Definition at line 49 of file DenseGslSolve.cpp.

49  { /*{{{*/
50 #ifdef _HAVE_GSL_
51  /*GSL Matrices and vectors: */
52  int s;
53  gsl_matrix_view a;
54  gsl_vector_view b,x;
55  gsl_permutation *p = NULL;
56 // for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
57 // for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
58  /*A will be modified by LU decomposition. Use copy*/
59  double* Acopy = xNew<double>(n*n);
60  xMemCpy(Acopy,A,n*n);
61 
62  /*Initialize gsl matrices and vectors: */
63  a = gsl_matrix_view_array (Acopy,n,n);
64  b = gsl_vector_view_array (B,n);
65  x = gsl_vector_view_array (X,n);
66 
67  /*Run LU and solve: */
68  p = gsl_permutation_alloc (n);
69  gsl_linalg_LU_decomp (&a.matrix, p, &s);
70  gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
71 
72  /*Clean up and assign output pointer*/
73  xDelete(Acopy);
74  gsl_permutation_free(p);
75 #endif
76 }
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
SolverxSeq
void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B, int n)
Definition: DenseGslSolve.cpp:49
xMemCpy
T * xMemCpy(T *dest, const T *src, unsigned int size)
Definition: MemOps.h:152
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38