DenseGslSolve.cpp File Reference

: solve dense matrix system with GSL library More...

#include "../../shared/shared.h"
#include "../../classes/Params/GenericParam.h"
#include "../../classes/Params/Parameters.h"
#include "../adolc/adolcincludes.h"
#include "../codipack/codipackincludes.h"
#include "./gslincludes.h"

Go to the source code of this file.

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

Detailed Description

: solve dense matrix system with GSL library

Definition in file DenseGslSolve.cpp.

Function Documentation

DenseGslSolve() [1/2]

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

Definition at line 25 of file DenseGslSolve.cpp.

                                                                         { /*{{{*/
 
   /*Intermediary: */
   IssmPDouble *X  = xNew<IssmPDouble>(n);
   SolverxSeq(X,A,B,n);
 
   /*allocate output pointers: */
   *pX=X;
}

DenseGslSolve() [2/2]

void DenseGslSolve	(	IssmPDouble **	px,
		IssmPDouble *	Kff,
		int	Kff_M,
		int	Kff_N,
		IssmPDouble *	pf,
		int	pf_M,
		Parameters *	parameters
	)

Definition at line 35 of file DenseGslSolve.cpp.

                                                                                                                         { /*{{{*/
 
   /*Intermediary: */
 
   if(Kff_N!=pf_M)_error_("Right hand side vector of size " << pf_M << ", when matrix is of size " << Kff_M << "-" << Kff_N << " !");
   if(Kff_M!=Kff_N)_error_("Stiffness matrix should be square!");
 
   IssmPDouble *x  = xNew<IssmPDouble>(Kff_N);
   SolverxSeq(x,Kff,pf,Kff_N);
 
   /*allocate output pointers: */
   *px=x;
}

SolverxSeq()

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

Definition at line 49 of file DenseGslSolve.cpp.

                                                                     { /*{{{*/
#ifdef _HAVE_GSL_
   /*GSL Matrices and vectors: */
   int              s;
   gsl_matrix_view  a;
   gsl_vector_view  b,x;
   gsl_permutation *p = NULL;
// for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
// for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
   /*A will be modified by LU decomposition. Use copy*/
   double* Acopy = xNew<double>(n*n);
   xMemCpy(Acopy,A,n*n);
 
   /*Initialize gsl matrices and vectors: */
   a = gsl_matrix_view_array (Acopy,n,n);
   b = gsl_vector_view_array (B,n);
   x = gsl_vector_view_array (X,n);
 
   /*Run LU and solve: */
   p = gsl_permutation_alloc (n);
   gsl_linalg_LU_decomp (&a.matrix, p, &s);
   gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
 
   /*Clean up and assign output pointer*/
   xDelete(Acopy);
   gsl_permutation_free(p);
#endif
}