#include <snowpack/snowpackCore/Solver.h>
#include <meteoio/MeteoIO.h>
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <math.h>

Classes
struct	pBLOCK

Macros
#define	GD_MALLOC(pointer, TYPE, N, MSG)

#define	GD_REALLOC(pointer, TYPE, N, MSG)

#define	GD_FREE(pointer)

#define	VD_AXPY_JUMP(N_B, N, JUMP, A, X, Y)

#define	VD_AXPY(N, A, X, Y)

#define	VD_AXPY_POS(N_B, N, POS, A, X, Y)

#define	VD_DOT_POS(N_B, N, POS, X, Y, RESULT)

#define	SD_CHUNK_REALLOC_SIZE 250

#define	SD_ALLOC_CHUNK(CHUNK, SIZE)

#define	SD_DESTROY_CHUNK(CHUNK)

#define	SD_MARKED (1<<30) /* An flag bit used for permutation. Use:(1<<15) on PC */

#define	SD_P_FIRST_COL_BLOCK(pMAT, pROW) (pMAT->pFirstColBlock + pROW->iColBlock)

#define	SD_P_SIZE_COL_BLOCK(pMAT, pROW) (pMAT->pSizeColBlock + pROW->iColBlock)

#define	SD_SEARCH_BLOCK_ROW(ROW, pROW_LOW, pROW_HIGH, pROW)

#define	SD_ROW( NUM, pMAT) ( (pMAT)->pRow[NUM] )

#define	SD_COL(pCOL) ( (pCOL)->Col )

#define	SD_N_ALLOC_COL 500

#define	SD_ALLOC_COL(N_COL, pMAT)

#define	SD_GET_COL(pCOL, pMAT)

#define	SD_FREE_COL_BLOCK_0(pCOL_BLOCK, pMAT) { (pCOL_BLOCK)->Next = (pMAT)->FreeColBlock; (pMAT)->FreeColBlock = (pCOL_BLOCK); }

#define	SD_FREE_COL_BLOCK(pCOL_BLOCK, pMAT) { SD_FREE_COL_BLOCK_0(pCOL_BLOCK, pMAT); pMAT->nColBlock--; }

#define	SD_ALLOC_COL_BLOCK(N_COL_BLOCK, pMAT)

#define	SD_N_ALLOC_COL_BLOCK 500

#define	SD_GET_COL_BLOCK(pCOL_BLOCK, pMAT)

#define	SD_FIND_COL(pROOT_COL, COL, ppCOL, FOUND)
	The SD_FIND_COL macro, accept pROOT_COL as the first node of the list to start the search of the node with value COL. If the node is found the variable FOUND is set to TRUE and ppCOL will point to this node, if not found FOUND is set to FALSE and ppCOL will point to the entry point. NOTE. For 2D meshes is better to not enable SPARSE_BINARY_TREE and use a linear list of column coefficients instead of a binary tree of column coefficients. More...

#define	SD_INSERT_COL(ppCOL, pCOL, COL)

#define	ERROR_SOLVER(MSG) { printf("++++Errror:%s:%d:%s\n", __FILE__, __LINE__, MSG); return(1); }

#define	USER_ERROR(MSG) { printf("++++Errror:gs_SolveMatrix:%s\n", MSG); return 1; }

#define	EXIT(MSG) { printf("++++Exit:%s:%d:%s\n", __FILE__, __LINE__, MSG); exit(1); }

#define	BLOCK_INIT(BLOCK, pCOL0, pSIZE) { BLOCK.pC0 = pCOL0; BLOCK.pSize = pSIZE; }

#define	BLOCK_NEXT(BLOCK) ( BLOCK.pC0++, BLOCK.pSize++ )

#define	BLOCK_C0(BLOCK) BLOCK.pC0[0]

#define	BLOCK_C1(BLOCK) (BLOCK.pC0[0]+BLOCK.pSize[0])

#define	BLOCK_SIZE(BLOCK) BLOCK.pSize[0]

#define	pC0_FIRST_COL( pROW) (pMatFirstColBlock + pROW->iColBlock)

#define	pSIZE_FIRST_COL(pROW) (pMatSizeColBlock + pROW->iColBlock)

#define	FIND_COL_BLOCK(pFIRST_BLK, COL, ppBLK, FOUND)

#define	FACT_SYM_MAT(MAT, N_ROW, N_COL)

#define	FIRST_BLOCK_ROW(pMAT) ( pMAT->pRowBlock )

#define	LAST_BLOCK_ROW(pMAT) ( pMAT->pRowBlock + pMAT->nRowBlock - 1 )

#define	SEARCH_ROW(ROW, pROW_LOW, pROW_HIGH, pROW)

#define	BLOCK_JUMP(nCOL0, pCOL0, pSIZE0, pCOL1, pSIZE1, JUMP)

#define	DIAGONAL(DIM, K) ( (K)(DIM) -( (K)((K)-1) )/2 )

#define	SEARCH_COL(COL, ROW, pMAT, pROW, FOUND, OFFSET)

#define	pFIRST_COL_BLOCK(pROW_BLOCK) ( ( (SD_TMP_ROW_BLOCK_DATA *)pROW_BLOCK)->Data.ColBlock )

Functions
int	Permute (int N, int Perm, double Vector)
	This function permute a vector, for a given permutation vector and compute the inverse permutation vector. Of course this is a trivial task if we have a second vector to store the new permuted vector, but we do not want to allocated extra memory so that we use a somewhat expensive algorithm which do not yet require a second storage vector. Here we assume that the dimension of permutation array is less than 2**31, because we use the 31th bit of each index of the permutation vector to store a flag, which tells us that we have already permuted that element. The algorithm is very simple, we take one element, look if we have already permuted that element, and if not we permute all elements connected to this element by a permutation cycle i.e. we shift right by one all elements of that cycle. ATTENTION: This function not only permute the given vector, but also compute the inverse permutation vector which is stored at the place of the permutation vector. More...

int	InvertMatrix (SD_BLOCK_MATRIX_DATA *pMat)
	This function is the kernel of the solution algorithm, and compute the LU triangular factorization of a block matrix where the blocks are so defined that there is no more fill-in during the factorizatio process. The LU factorization of the matrix substitute the original matrix data. The LU factorization of a matrix is very similar to the inverse of the matrix and from this fact we have derived the name for this function. More...

int	InverseMatrixVector (SD_BLOCK_MATRIX_DATA pMat, double X)
	This function perform the backward and forward substitution for a LU factorized block matrix. Of course this operation is the same as multipying the inverse matrix with a vector and inf fact the implementation does not much differs from a matrix vector multiplication. More...

int	BuildSparseConFormat (SD_CON_MATRIX_DATA pMat, int pRowStart0, int *pColumn0)
	This function build a compressed sparse matrix format which is the common used one. The value pRowStart[Row] specifies the offset in the array pColumn[] where the column coefficients for the row: Row are starting to be defined. The column coefficients for a row are always sorted. Of course the array pRowStart has a dimension of plus one with respect to the matrix dimension. ATTENTION: Here we use a FORTRAN notation, all values in the arrays pRowStart and pColumn are incremented by one. More...

static void	RunMmd (int neqns, int xadj, int adjncy, int invp, int perm, int delta, int head, int qsize, int nsize, int list, int marker, int maxint, int ncsub)
	RunMmd – multiple minimum external degree purpose – this routine implements the minimum degree algorithm. it makes use of the implicit representation of elimination graphs by quotient graphs, and the notion of indistinguishable nodes. It also implements the modifications by multiple elimination and minimum external degree. Caution – the adjacency vector adjncy will be destroyed. Input parameters –. More...

int	ComputeTmpConMatrix (SD_CON_MATRIX_DATA pMat0, SD_TMP_CON_MATRIX_DATA pMat)
	This routine set up a temporary adjacency block matrix format for the purpose to speed-up the symbolic factorization process and is called after the permutation vector has been computed. From this step all indices for rows and columns are permuted indices. The row blocks are a result of the mmd algorithm and so need not to be computed. Here we have to compute the column blocks for the upper matrix including the diagonal. Of course the matrix is still sparse because at this time the fill-in is still unknown but it is however possible that blocks with size larger than one will appear. NOTE: All data allocated in pMat0 is released. More...

int	ComputeFillIn (SD_TMP_CON_MATRIX_DATA *pMat)
	This function compute the fill-in by factorizing symbolically the matrix. More...

int	ComputeBlockMatrix (SD_TMP_CON_MATRIX_DATA pTmpMat, SD_BLOCK_MATRIX_DATA pMat)
	Here, we define the matrix as used by the numerical factorization algorithm. If the multiplicity factor is not 1 at this time we define the true row and column blocks. More...

int	AllocateConData (int Dim, SD_CON_MATRIX_DATA *pMat)
	Allocate the matrix row data in order to store the connectivity matrix data. More...

int	ReleaseConMatrix (SD_CON_MATRIX_DATA *pMat)

int	ds_Initialize (const int &MatDim, SD_MATRIX_DATA **ppMat)
	This is the first function to be called prior to begin to work with any matrix. The user must gives the number of equations i.e. the dimension of the matrix, the type of matrix i.e. a symmetric or unsymmetric matrix and the multiplicity. The multiplicity specifies that each coefficient of the matrix is actually a square matrix and each equation and unknown is actually a set of equations and unknowns all of dimension equal to the defined multiplicity factor. Thus the true dimension of the linear system is given by the specified matrix dimension times the multiplicity. Of course a multiplicity value of one is the most general case, but sometines especially by vector field computations the multiplicity is simply given by the dimension of the vector to be computed and this allow to speed up many integer operations. Only define a multiplicity factor greather than one, if the vector field components, or a subset, are fully coupled to eachother. In fact, by definition of a multiplicity we always reserve memory for the full coupled system among each vector component, thus the memory requirement increase quadratically with the defined multiplicity factor. The function return a pointer to an opaque data type as matrix identifier. More...

void	FACT_SYM_MAT_BLOCK (int N_PIVOT, int TOT_ROW, int N_ROW, int N_COL, double MAT0, int DIM0, double MAT1, int DIM1, int N_BLOCK, int N, int JUMP)

int	ComputePermutation (SD_CON_MATRIX_DATA *pMat)
	This function compute the permutation array and thus run the mmd algorithm. More...

int	SymbolicFact (SD_MATRIX_DATA *pMat)

bool	ds_Solve (const SD_MATRIX_WHAT &Code, SD_MATRIX_DATA pMat, double X)
	This function calls the solver itself. More...

int	ds_AssembleMatrix (SD_MATRIX_DATA pMat0, const int &nEq, int Eq[], const int &Dim, const double ElMat)
	This function assemble the element square matrix [ElMat] for one element with nEqM x nEqM real coefficients in to the global matrix [A]. It must be called for each (finite) element after the element connectivity have been assembled and the matrix symbolic factorized. To perform this task we also require the element incidences. The variable: Dim specifies the dimension of the matrix: Mat which is not required to be equal to the numer of element incidences: nEq. More...

static void	MmdUpdate (int ehead, int neqns, int xadj, int adjncy, int delta, int mdeg, int head, int forward, int backward, int qsize, int list, int marker, int maxint, int tag)
	MmdUpdate -— multiple minimum degree update purpose – this routine updates the degrees of nodes after a multiple elimination step. input parameters –. More...

static void	MmdElimin (int mdeg_node, int xadj, int adjncy, int head, int forward, int backward, int qsize, int list, int marker, int maxint, int tag)
	MmdElimin – multiple minimum degree elimination Purpose – This routine eliminates the node mdeg_node of minimum degree from the adjacency structure, which is stored in the quotient graph format. It also transforms the quotient graph representation of the elimination graph. Input parameters –. More...

static int	MmdInit (int neqns, int xadj, int head, int forward, int backward, int qsize, int list, int *marker)
	MmdInit – mult minimum degree initialization purpose – this routine performs initialization for the multiple elimination version of the minimum degree algorithm. input parameters –. More...

static void	MmdNumbering (int neqns, int perm, int invp, int qsize, int nsize)
	MmdNumbering – multi minimum degree numbering purpose – this routine performs the final step in producing the permutation and inverse permutation vectors in the multiple elimination version of the minimum degree ordering algorithm. input parameters –. More...

void	MERGE_COL_BLOCK (SD_COL_BLOCK_DATA pCOL0, SD_COL_BLOCK_DATA ppCOL1, SD_TMP_CON_MATRIX_DATA pMAT)

int	ReleaseBlockMatrix (SD_BLOCK_MATRIX_DATA *pMat)
	Release all the data allocated for the numerical factorization algorithm. More...

int	ds_DefineConnectivity (SD_MATRIX_DATA *pMat0, const int &nEq, int Eq[], const int &nEl, const int &Dim)
	This function assemble the element connnectivity for one or more elements in order to build a sparse matrix format. Of course we only store the upper part of the connectivity matrix because we only consider structure symmetric matrices. More...

Variables
static bool	gd_MemErr

Macro Definition Documentation

◆ BLOCK_C0

#define BLOCK_C0 ( BLOCK ) BLOCK.pC0[0]

◆ BLOCK_C1

#define BLOCK_C1 ( BLOCK ) (BLOCK.pC0[0]+BLOCK.pSize[0])

◆ BLOCK_INIT

#define BLOCK_INIT	(	BLOCK,
		pCOL0,
		pSIZE
	)	{ BLOCK.pC0 = pCOL0; BLOCK.pSize = pSIZE; }

◆ BLOCK_JUMP

#define BLOCK_JUMP	(	nCOL0,
		pCOL0,
		pSIZE0,
		pCOL1,
		pSIZE1,
		JUMP
	)

Value:

{  int i_, *pCol0_, *pCol1_, *pSize0_,  *pSize1_, Size_, Col1_0_, Col1_1_;                     \
   pCol0_  = pCOL0;                                                                            \
   pCol1_  = pCOL1;                                                                            \
   pSize0_ = pSIZE0;                                                                           \
   pSize1_ = pSIZE1;                                                                           \
   Col1_0_ = pCol1_[0];                                                                        \
   Col1_1_ = Col1_0_ + pSize1_[0];                                                             \
   for(i_=0; i_<nCOL0; i_++)                                                                   \
   {  Size_   = 0;                                                                             \
      while( pCol1_[0] + pSize1_[0] < pCol0_[0] )                                              \
      {  Size_ += Col1_1_ - Col1_0_;                                                           \
         Col1_1_  = ( Col1_0_ = (++pCol1_)[0] ) + (++pSize1_)[0];   }                          \
      JUMP[i_] = Size_ + pCol0_[0] - Col1_0_;                                                  \
      Col1_0_  = (pCol0_++)[0] + (pSize0_++)[0];                                               \
   }                                                                                           \
}

◆ BLOCK_NEXT

#define BLOCK_NEXT ( BLOCK ) ( BLOCK.pC0++, BLOCK.pSize++ )

◆ BLOCK_SIZE

#define BLOCK_SIZE ( BLOCK ) BLOCK.pSize[0]

◆ DIAGONAL

#define DIAGONAL	(	DIM,
		K
	)	( (K)(DIM) -( (K)((K)-1) )/2 )

◆ ERROR_SOLVER

#define ERROR_SOLVER ( MSG ) { printf("++++Errror:%s:%d:%s\n", __FILE__, __LINE__, MSG); return(1); }

◆ EXIT

#define EXIT ( MSG ) { printf("++++Exit:%s:%d:%s\n", __FILE__, __LINE__, MSG); exit(1); }

◆ FACT_SYM_MAT

#define FACT_SYM_MAT	(	MAT,
		N_ROW,
		N_COL
	)

Value:

{                                                                                              \
   if ( N_ROW>1 ) {                                                                            \
   const int m_n_1=N_COL-N_ROW+1;                                                              \
   double *Mat_k=MAT;                                                                           \
   for ( int n_k=N_COL; n_k>=m_n_1; n_k-- )                                                    \
   {  double Pivot = 1./(*Mat_k);                                                               \
      double *Mat_i = Mat_k++;                                                                  \
      for ( int n_i=n_k; n_i>m_n_1; Mat_k++ )                                                  \
      {  Mat_i += n_i--;  VD_AXPY(n_i, -(*Mat_k)*Pivot, Mat_k, Mat_i);  }                      \
      Mat_k += m_n_1 - 1;                                                                      \
   }                                                                                           \
  }                                                                                            \
}

◆ FIND_COL_BLOCK

#define FIND_COL_BLOCK	(	pFIRST_BLK,
		COL,
		ppBLK,
		FOUND
	)

Value:

{                                                                                              \
   FOUND      = 0;                                                                             \
   SD_COL_BLOCK_DATA *pB_ = (pFIRST_BLK);                                                      \
   int Col0_  = COL+1;                                                                      \
   int Col1_  = COL-1;                                                                      \
   while ( pB_ )                                                                               \
   {  if      ( Col1_ >  pB_->Col1  ) { ppBLK = &pB_->Next; pB_ = pB_->Next; }                 \
      else if ( Col0_ >= pB_->Col0  )                                                          \
      {  FOUND = 1;                                                                            \
         if      ( Col0_==pB_->Col0 ) pB_->Col0 = COL;                                         \
         else if ( Col1_==pB_->Col1 ) pB_->Col1 = COL;                                         \
         break;                                                                                \
      }                                                                                        \
      else   break;                                                                            \
   }                                                                                           \
}

◆ FIRST_BLOCK_ROW

#define FIRST_BLOCK_ROW ( pMAT ) ( pMAT->pRowBlock )

◆ GD_FREE

#define GD_FREE ( pointer )

Value:

{                                                                                              \
        if ( pointer ) {                                                                       \
                free ( (char*) pointer );                                                      \
                pointer = NULL;                                                                \
        }                                                                                      \
}                                                                                              \

◆ GD_MALLOC

#define GD_MALLOC	(	pointer,
		TYPE,
		N,
		MSG
	)

Value:

{                                                                                              \
        pointer = (TYPE *)malloc( sizeof(TYPE)*(N+1) );                                        \
        if ( pointer  ) {                                                                      \
                gd_MemErr = false;                                                             \
        } else {                                                                               \
                gd_MemErr = true; fprintf(stderr, "\n+++++ %s: %s\n",  "NO SPACE TO ALLOCATE", MSG);               \
        }                                                                                      \
}

◆ GD_REALLOC

#define GD_REALLOC	(	pointer,
		TYPE,
		N,
		MSG
	)

Value:

{                                                                                              \
        if ( pointer )  {                                                                      \
                pointer = (TYPE *)realloc( (char*)pointer, sizeof(TYPE)*(N+1) );               \
                if ( pointer  ) {                                                                      \
                        gd_MemErr = false;                                                             \
                } else {                                                                               \
                        gd_MemErr = true; fprintf(stderr, "\n+++++ %s: %s\n", "NO SPACE TO REALLOCATE", MSG);               \
                }                                                                                   \
        } else {                                                                               \
                GD_MALLOC(  pointer, TYPE, N, MSG );                                           \
        }                                                                                      \
}

◆ LAST_BLOCK_ROW

#define LAST_BLOCK_ROW ( pMAT ) ( pMAT->pRowBlock + pMAT->nRowBlock - 1 )

◆ pC0_FIRST_COL

#define pC0_FIRST_COL ( pROW ) (pMatFirstColBlock + pROW->iColBlock)

◆ pFIRST_COL_BLOCK

#define pFIRST_COL_BLOCK ( pROW_BLOCK ) ( ( (SD_TMP_ROW_BLOCK_DATA *)pROW_BLOCK)->Data.ColBlock )

◆ pSIZE_FIRST_COL

#define pSIZE_FIRST_COL ( pROW ) (pMatSizeColBlock + pROW->iColBlock)

◆ SD_ALLOC_CHUNK

#define SD_ALLOC_CHUNK	(	CHUNK,
		SIZE
	)

Value:

{  if ( CHUNK.nChunks >= CHUNK.pChunksSize )                                                   \
   {  CHUNK.pChunksSize += SD_CHUNK_REALLOC_SIZE;                                              \
      GD_REALLOC( CHUNK.pChunks, char*, CHUNK.pChunksSize, "Chunk pointer Data" );             \
   }                                                                                           \
   GD_MALLOC( CHUNK.pChunks[ CHUNK.nChunks ], char, SIZE, "Chunk Data" );                      \
   CHUNK.TotChunkSize += SIZE;                                                                 \
   CHUNK.nChunks++;                                                                            \
}

◆ SD_ALLOC_COL

#define SD_ALLOC_COL	(	N_COL,
		pMAT
	)

Value:

{  SD_ALLOC_CHUNK((pMAT)->PoolCol, sizeof(SD_COL_DATA)*N_COL);                                 \
   (pMAT)->FreeCol = ( SD_COL_DATA * ) (pMAT)->PoolCol.pChunks[ (pMAT)->PoolCol.nChunks-1 ];   \
   (pMAT)->nFreeCol = N_COL;                                                                   \
}

◆ SD_ALLOC_COL_BLOCK

#define SD_ALLOC_COL_BLOCK	(	N_COL_BLOCK,
		pMAT
	)

Value:

{  SD_COL_BLOCK_DATA *pColBlock_;                                                              \
   SD_ALLOC_CHUNK((pMAT)->PoolColBlock, sizeof(SD_COL_BLOCK_DATA)*N_COL_BLOCK);                \
   pColBlock_ = ( SD_COL_BLOCK_DATA * )                                                        \
                (pMAT)->PoolColBlock.pChunks[ (pMAT)->PoolColBlock.nChunks-1 ];                \
   for(int i_=N_COL_BLOCK; 0<i_--; pColBlock_++) SD_FREE_COL_BLOCK_0(pColBlock_, pMAT);        \
}

◆ SD_CHUNK_REALLOC_SIZE

#define SD_CHUNK_REALLOC_SIZE 250

◆ SD_COL

#define SD_COL ( pCOL ) ( (pCOL)->Col )

◆ SD_DESTROY_CHUNK

#define SD_DESTROY_CHUNK ( CHUNK )

Value:

{  int i_;                                                                                     \
   for(i_=0; i_<CHUNK.nChunks; i_++) GD_FREE(CHUNK.pChunks[i_]);                               \
   GD_FREE(CHUNK.pChunks);                                                                     \
   CHUNK.TotChunkSize = 0;                                                                     \
}

◆ SD_FIND_COL

#define SD_FIND_COL	(	pROOT_COL,
		COL,
		ppCOL,
		FOUND
	)

Value:

{  SD_COL_DATA *pC_ ;                                                                          \
   FOUND   = 0;                                                                                \
   pC_     = (pROOT_COL);                                                                      \
   ppCOL   = &(pROOT_COL);                                                                     \
   while ( pC_ )                                                                               \
   {  if ( COL > SD_COL(pC_)  )  { ppCOL = &pC_->Next; pC_ = pC_->Next;  }                     \
      else { if ( COL == SD_COL(pC_) ) FOUND = 1;  break;  }                                   \
   }                                                                                           \
}

The SD_FIND_COL macro, accept pROOT_COL as the first node of the list to start the search of the node with value COL. If the node is found the variable FOUND is set to TRUE and ppCOL will point to this node, if not found FOUND is set to FALSE and ppCOL will point to the entry point. NOTE. For 2D meshes is better to not enable SPARSE_BINARY_TREE and use a linear list of column coefficients instead of a binary tree of column coefficients.

◆ SD_FREE_COL_BLOCK

#define SD_FREE_COL_BLOCK	(	pCOL_BLOCK,
		pMAT
	)	{ SD_FREE_COL_BLOCK_0(pCOL_BLOCK, pMAT); pMAT->nColBlock--; }

◆ SD_FREE_COL_BLOCK_0

#define SD_FREE_COL_BLOCK_0	(	pCOL_BLOCK,
		pMAT
	)	{ (pCOL_BLOCK)->Next = (pMAT)->FreeColBlock; (pMAT)->FreeColBlock = (pCOL_BLOCK); }

◆ SD_GET_COL

#define SD_GET_COL	(	pCOL,
		pMAT
	)

Value:

{  if  ( !(pMAT)->nFreeCol )  SD_ALLOC_COL(SD_N_ALLOC_COL, pMAT);                              \
   pCOL = (pMAT)->FreeCol++; (pMAT)->nFreeCol--;                                               \
}

◆ SD_GET_COL_BLOCK

#define SD_GET_COL_BLOCK	(	pCOL_BLOCK,
		pMAT
	)

Value:

{  if  ( !(pMAT)->FreeColBlock )  SD_ALLOC_COL_BLOCK(SD_N_ALLOC_COL_BLOCK, pMAT);              \
   pCOL_BLOCK = (pMAT)->FreeColBlock; (pMAT)->FreeColBlock = (pMAT)->FreeColBlock->Next;       \
   pMAT->nColBlock++;                                                                          \
}

◆ SD_INSERT_COL

#define SD_INSERT_COL	(	ppCOL,
		pCOL,
		COL
	)

Value:

{  pCOL->Col  = COL;                                                                           \
   pCOL->Next = *ppCOL;                                                                        \
   *ppCOL     = pCOL;                                                                          \
}

◆ SD_MARKED

#define SD_MARKED (1<<30) /* An flag bit used for permutation. Use:(1<<15) on PC */

◆ SD_N_ALLOC_COL

#define SD_N_ALLOC_COL 500

◆ SD_N_ALLOC_COL_BLOCK

#define SD_N_ALLOC_COL_BLOCK 500

◆ SD_P_FIRST_COL_BLOCK

#define SD_P_FIRST_COL_BLOCK	(	pMAT,
		pROW
	)	(pMAT->pFirstColBlock + pROW->iColBlock)

◆ SD_P_SIZE_COL_BLOCK

#define SD_P_SIZE_COL_BLOCK	(	pMAT,
		pROW
	)	(pMAT->pSizeColBlock + pROW->iColBlock)

◆ SD_ROW

#define SD_ROW	(	NUM,
		pMAT
	)	( (pMAT)->pRow[NUM] )

◆ SD_SEARCH_BLOCK_ROW

#define SD_SEARCH_BLOCK_ROW	(	ROW,
		pROW_LOW,
		pROW_HIGH,
		pROW
	)

Value:

{  SD_ROW_BLOCK_DATA *low_, *high_, *mid_;                                                     \
   low_ = pROW_LOW; high_ = pROW_HIGH;                                                         \
   while( low_<=high_ )                                                                        \
   {  mid_ = low_ + ( high_ - low_ ) / 2;                                                      \
      if      ( ROW < mid_->Row0 ) high_ = mid_ - 1;                                           \
      else if ( ROW > mid_->Row1 ) low_  = mid_ + 1;                                           \
      else { pROW=mid_;  break;  }                                                             \
   }                                                                                           \
}

◆ SEARCH_COL

#define SEARCH_COL	(	COL,
		ROW,
		pMAT,
		pROW,
		FOUND,
		OFFSET
	)

Value:

{  int *col_, *size_;                                                   \
   col_     = SD_P_FIRST_COL_BLOCK(pMAT,pROW);                                                 \
   size_    = SD_P_SIZE_COL_BLOCK( pMAT,pROW);                                                 \
   const int delta_   = ROW - pROW->Row0;                                                                \
   OFFSET   = pROW->iFloat + DIAGONAL(pROW->nCol, delta_);                                     \
   {  ++col_;                                                                                  \
      for(int i_=pROW->nColBlock-1; (i_--)>0; OFFSET += size_[0], col_++, size_++)                 \
      {  if ( COL < col_[0] )  break;  }                                                       \
       --col_;                                                                                 \
      if ( COL >= col_[0]+size_[0] ) {  FOUND = 0; }                                           \
      else                           {  FOUND = 1;  OFFSET += COL - col_[0] - delta_;  }       \
   }                                                                                           \
}

◆ SEARCH_ROW

#define SEARCH_ROW	(	ROW,
		pROW_LOW,
		pROW_HIGH,
		pROW
	)

Value:

{  SD_ROW_BLOCK_DATA *low_, *high_, *mid_;                                                     \
   low_ = pROW_LOW; high_ = pROW_HIGH;                                                         \
   while( low_<=high_ )                                                                        \
   {  mid_ = low_ + ( high_ - low_ ) / 2;                                                      \
      if      ( ROW < mid_->Row0 ) high_ = mid_ - 1;                                           \
      else if ( ROW > mid_->Row1 ) low_  = mid_ + 1;                                           \
      else { pROW=mid_;  break;  }                                                             \
   }                                                                                           \
}

◆ USER_ERROR

#define USER_ERROR ( MSG ) { printf("++++Errror:gs_SolveMatrix:%s\n", MSG); return 1; }

◆ VD_AXPY

#define VD_AXPY	(	N,
		A,
		X,
		Y
	)

Value:

/* Y[] = A*X[] + Y[] */                                      \
{  double  a_, *x_, *y_;                                                               \
   int  k_;                                                                           \
   for (x_=X, y_=Y, a_=A, k_=N; 0<k_--; ) *y_++ += (a_)*(*x_++) ;                              \
}

◆ VD_AXPY_JUMP

#define VD_AXPY_JUMP	(	N_B,
		N,
		JUMP,
		A,
		X,
		Y
	)

Value:

/* Y[] = A*X[] + Y[] BLOCK-WISE IN Y */           \
{  double  a_=A, *x_=X, *y_=Y+JUMP[0];                                                         \
   int  n_;                                                                                    \
   for (n_=0; n_<N_B; y_+= JUMP[++n_])                                                         \
   for (int k_=N[n_]; 0<k_--;  ) *y_++ += (a_)*(*x_++) ;                                       \
}

◆ VD_AXPY_POS

#define VD_AXPY_POS	(	N_B,
		N,
		POS,
		A,
		X,
		Y
	)

Value:

/* Y[] = A*X[] + Y[] BLOCK-WISE IN Y */             \
{  double  a_, *x_, *y_;                                                               \
   int  n_, k_;                                                                       \
   for (x_=X, y_=Y+POS[0], a_=A, n_=0; n_<N_B; y_= Y+POS[++n_])                                \
   for (k_=N[n_]; 0<k_--;  ) *y_++ += (a_)*(*x_++) ;                                           \
}

◆ VD_DOT_POS

#define VD_DOT_POS	(	N_B,
		N,
		POS,
		X,
		Y,
		RESULT
	)

Value:

/* RESULT = X[]*Y[]  BLOCK-WISE IN Y */          \
{  double  r_, *x_, *y_;                                                               \
   int  n_, k_;                                                                       \
   for (x_=X, y_=Y+POS[0], r_=0.0, n_=0; n_<N_B; y_= Y+POS[++n_])                              \
   for (k_=N[n_]; 0<k_--;  ) r_ += (*x_++)*(*y_++) ;                                           \
   RESULT = r_;                                                                                \
}

Function Documentation

◆ AllocateConData()

int AllocateConData	(	int	Dim,
		SD_CON_MATRIX_DATA *	pMat
	)

inline

Allocate the matrix row data in order to store the connectivity matrix data.

Parameters

Dim	int
pMat	SD_CON_MATRIX_DATA

Returns: int

◆ BuildSparseConFormat()

int BuildSparseConFormat	(	SD_CON_MATRIX_DATA *	pMat,
		int *	pRowStart0,
		int *	pColumn0
	)

This function build a compressed sparse matrix format which is the common used one. The value pRowStart[Row] specifies the offset in the array pColumn[] where the column coefficients for the row: Row are starting to be defined. The column coefficients for a row are always sorted. Of course the array pRowStart has a dimension of plus one with respect to the matrix dimension. ATTENTION: Here we use a FORTRAN notation, all values in the arrays pRowStart and pColumn are incremented by one.

Parameters

pMat	SD_CON_MATRIX_DATA
pRowStart0	int
pColumn0	int

Returns: int

◆ ComputeBlockMatrix()

int ComputeBlockMatrix	(	SD_TMP_CON_MATRIX_DATA *	pTmpMat,
		SD_BLOCK_MATRIX_DATA *	pMat
	)

Here, we define the matrix as used by the numerical factorization algorithm. If the multiplicity factor is not 1 at this time we define the true row and column blocks.

Parameters

pTmpMat	SD_TMP_CON_MATRIX_DATA
pMat	SD_BLOCK_MATRIX_DATA

Returns: int

◆ ComputeFillIn()

int ComputeFillIn ( SD_TMP_CON_MATRIX_DATA * pMat )

This function compute the fill-in by factorizing symbolically the matrix.

Parameters

pMat	SD_TMP_CON_MATRIX_DATA

Returns: int

◆ ComputePermutation()

int ComputePermutation ( SD_CON_MATRIX_DATA * pMat )

inline

This function compute the permutation array and thus run the mmd algorithm.

Parameters

pMat	SD_CON_MATRIX_DATA

Returns: int

◆ ComputeTmpConMatrix()

int ComputeTmpConMatrix	(	SD_CON_MATRIX_DATA *	pMat0,
		SD_TMP_CON_MATRIX_DATA *	pMat
	)

This routine set up a temporary adjacency block matrix format for the purpose to speed-up the symbolic factorization process and is called after the permutation vector has been computed. From this step all indices for rows and columns are permuted indices. The row blocks are a result of the mmd algorithm and so need not to be computed. Here we have to compute the column blocks for the upper matrix including the diagonal. Of course the matrix is still sparse because at this time the fill-in is still unknown but it is however possible that blocks with size larger than one will appear. NOTE: All data allocated in pMat0 is released.

Parameters

pMat0	SD_CON_MATRIX_DATA
pMat	SD_TMP_CON_MATRIX_DATA

Returns: int

◆ ds_AssembleMatrix()

int ds_AssembleMatrix	(	SD_MATRIX_DATA *	pMat0,
		const int &	nEq,
		int	Eq[],
		const int &	Dim,
		const double *	ElMat
	)

This function assemble the element square matrix [ElMat] for one element with nEq*M x nEq*M real coefficients in to the global matrix [A]. It must be called for each (finite) element after the element connectivity have been assembled and the matrix symbolic factorized. To perform this task we also require the element incidences. The variable: Dim specifies the dimension of the matrix: Mat which is not required to be equal to the numer of element incidences: nEq.

If a multiplicity factor M greather than 1 has been defined the numerical values in the element matrix [ElMat] must be forall vector components clustered together. E.g. for a multiplicity of 3 i.e. a 3D vector field, the 3x3 left-upper submatrix of [ElMat] must represent the coupling between the 3 vector field components. To perform this task we also newly require the element connectivity. The list of element equations should be the same or a subset of the list used previously when the matrix connectivity has been defined with a call to: ds_DefineConnectivity(). ATTENTION: no error is detected if some of the the element connectivity defined when calling ds_AssembleLocalMatrix() are not included in those previously specified when calling ds_DefineConnectivity() (ie. if the specified incidences have not been previously defined).

If the matrix [A] has been declared as symmetric only the upper triangular part of [ElMat], i.e. only ElMat[i][j] = ElMat[Dim*i+j] with i<=j and i,j = 0...M*nEq-1 are used and need to be defined. In the unsymmetric case all M*nEq x M*nEq coefficients are used. It is assumed that in the calling program the array [ElMat] is dimensioned as ElMat[...][Dim] NOTICE: Of course the parameter Dim can be greater than M*nEq: Dim >= M*nEq After all element matrices have been assembled the matrix [A] = [L][U] can be factorised in to the lower [L] and upper [U] tringular matrices by calling ds_Solve(NumericFactorize, ).

Parameters

[in]	pMat0	pointer to the matrix [A] opaque data returned by ds_Initialize()
[in]	nEq	no. of equations for one element forming a crique
[in]	Eq	Element list of equations for one element.
[in]	Dim	first dimension of the 2D-array ElMat[][Dim]
[in]	ElMat	element square matrix to be assembled in the matrix [A]

◆ ds_DefineConnectivity()

int ds_DefineConnectivity	(	SD_MATRIX_DATA *const	pMat0,
		const int &	nEq,
		int	Eq[],
		const int &	nEl,
		const int &	Dim
	)

This function assemble the element connnectivity for one or more elements in order to build a sparse matrix format. Of course we only store the upper part of the connectivity matrix because we only consider structure symmetric matrices.

The next function represents the computational kernel of this direct solver. Its functionality depends on the input parameter: Code whose meanings are given below. Generally once the matrix as been defined and the element connectivity assembled, the user first perform a symbolic factorization process. After assembling the numerical matrix, a numerical factorization step follows together with a back- and for-ward substitution to compute a numerical solution of the linear system. The matrix data can be resetted in order to solve other linear systems having the same matrix connectivity without the need to newly symbolic factorize the matrix. NOTE: If a multiplicity factor has been defined we assume that the vector components are clustered together in the vector pVec.

it is needed for defining the system (matrix) connectivity i.e. the non-zero coefficients of the matrix [A] i.e. which equation is connected to which one. For each (finite) element we have to specifies a list of equations. Here, we assume that all equations in the list are connected to eachother and thus lead to non-zero coefficients in the matrix i.e. the element list of equations form a crique This function is generally called for a single element (nEl = 1) or for a group of nEl elements having criques of equal dimension. To contemporary define the element connectivity for more elements set nEl>1 and store the list equations for each element in the array: Eq[][Dim]. Of course nEl<=Dim must holds. If a multiplicity factor greather than 1 has been defined we have only to define the element connectivity for the first representative equations or first vector field component i.e. as we would do for a scalar field, and thus all equations in the list must have a number less than the specified matrix dimension E.g. for an element with 4 nodes and multiplicity 3 we have a total true number of 12 equations forming a crique but only 4 = ( 12 / 3 ) equations must be specified and the values must not exceed the specified dimension of the matrix [A]. After the list of equations for all elements have been specified the matrix [A] should be symbolically factorized by calling the function: ds_Solve(SymbolicFactorize, ... ). NOTE: Except the definition of the multiplicity in ds_Initialize(), all steps performed to define the structure of matrix [A] are stricktly independent from the multiplicity

Parameters

[in]	pMat0	pointer to the matrix [A] opaque data returned by ds_Initialize()
[in]	nEq	No. of equations for one element forming a crique
[in]	Eq	Element list of equations for more elements with equal no. of eqs.
[in]	nEl	No. of elements ( 0 <= i "<" nEq ; 0 <= e "<" nEl )
[in]	Dim	first dimension of the 2D-array Eq[][Dim]

◆ ds_Initialize()

int ds_Initialize	(	const int &	MatDim,
		SD_MATRIX_DATA **	ppMat
	)

This is the first function to be called prior to begin to work with any matrix. The user must gives the number of equations i.e. the dimension of the matrix, the type of matrix i.e. a symmetric or unsymmetric matrix and the multiplicity. The multiplicity specifies that each coefficient of the matrix is actually a square matrix and each equation and unknown is actually a set of equations and unknowns all of dimension equal to the defined multiplicity factor. Thus the true dimension of the linear system is given by the specified matrix dimension times the multiplicity. Of course a multiplicity value of one is the most general case, but sometines especially by vector field computations the multiplicity is simply given by the dimension of the vector to be computed and this allow to speed up many integer operations. Only define a multiplicity factor greather than one, if the vector field components, or a subset, are fully coupled to eachother. In fact, by definition of a multiplicity we always reserve memory for the full coupled system among each vector component, thus the memory requirement increase quadratically with the defined multiplicity factor. The function return a pointer to an opaque data type as matrix identifier.

Parameters

MatDim	dimension of the matrix [A]. The true number of equations and unknowns is given by: MatDim * Multiplicity
ppMat	A pointer to an opaque data type storing data related to the matrix [A]

◆ ds_Solve()

bool ds_Solve	(	const SD_MATRIX_WHAT &	Code,
		SD_MATRIX_DATA *	pMat,
		double *	pX
	)

This function calls the solver itself.

Parameters

[in]	Code	functionlity code defined above
[in]	pMat	pointer to the matrix [A] opaque data
[in]	pX	right hand side vector {B} to be overwritten by the solution vector {X}: B[i] := X[i]

Returns: false whenever the solve produced NaNs in the solution vector

◆ FACT_SYM_MAT_BLOCK()

void FACT_SYM_MAT_BLOCK	(	int	N_PIVOT,
		int	TOT_ROW,
		int	N_ROW,
		int	N_COL,
		double *	MAT0,
		int	DIM0,
		double *	MAT1,
		int	DIM1,
		int	N_BLOCK,
		int *	N,
		int *	JUMP
	)

inline

◆ InverseMatrixVector()

int InverseMatrixVector	(	SD_BLOCK_MATRIX_DATA *	pMat,
		double *	X
	)

This function perform the backward and forward substitution for a LU factorized block matrix. Of course this operation is the same as multipying the inverse matrix with a vector and inf fact the implementation does not much differs from a matrix vector multiplication.

Parameters

pMat	SD_BLOCK_MATRIX_DATA
X	double

Returns: int

◆ InvertMatrix()

int InvertMatrix ( SD_BLOCK_MATRIX_DATA * pMat )

This function is the kernel of the solution algorithm, and compute the LU triangular factorization of a block matrix where the blocks are so defined that there is no more fill-in during the factorizatio process. The LU factorization of the matrix substitute the original matrix data. The LU factorization of a matrix is very similar to the inverse of the matrix and from this fact we have derived the name for this function.

Parameters

pMat	SD_BLOCK_MATRIX_DATA

Returns: int

◆ MERGE_COL_BLOCK()

void MERGE_COL_BLOCK	(	SD_COL_BLOCK_DATA *	pCOL0,
		SD_COL_BLOCK_DATA **	ppCOL1,
		SD_TMP_CON_MATRIX_DATA *	pMAT
	)

inline

◆ MmdElimin()

static void MmdElimin	(	int	mdeg_node,
		int *	xadj,
		int *	adjncy,
		int *	head,
		int *	forward,
		int *	backward,
		int *	qsize,
		int *	list,
		int *	marker,
		int	maxint,
		int	tag
	)

static

MmdElimin – multiple minimum degree elimination Purpose – This routine eliminates the node mdeg_node of minimum degree from the adjacency structure, which is stored in the quotient graph format. It also transforms the quotient graph representation of the elimination graph. Input parameters –.

Parameters

mdeg_node	– node of minimum degree.
maxint	– estimate of maximum representable (short) integer.
tag	– tag value. Updated parameters –
(xadj,adjncy)	– updated adjacency structure.
(head,forward,backward)	– degree doubly linked structure.
qsize	– size of supernode.
marker	– marker vector.
list	– temporary linked list of eliminated nabors.

◆ MmdInit()

static int MmdInit	(	int	neqns,
		int *	xadj,
		int *	head,
		int *	forward,
		int *	backward,
		int *	qsize,
		int *	list,
		int *	marker
	)

static

MmdInit – mult minimum degree initialization purpose – this routine performs initialization for the multiple elimination version of the minimum degree algorithm. input parameters –.

Parameters

neqns	– number of equations.
xadj	– adjacency structure. output parameters –
(head,forward,backward)	– degree doubly linked structure.
qsize	– size of supernode ( initialized to one).
list	– linked list.
marker	– marker vector.

◆ MmdNumbering()

static void MmdNumbering	(	int	neqns,
		int *	perm,
		int *	invp,
		int *	qsize,
		int *	nsize
	)

static

MmdNumbering – multi minimum degree numbering purpose – this routine performs the final step in producing the permutation and inverse permutation vectors in the multiple elimination version of the minimum degree ordering algorithm. input parameters –.

Parameters

neqns	– number of equations.
qsize	– size of supernodes at elimination. updated parameters –
invp	– inverse permutation vector. on input, if qsize[node] = 0, then node has been output parameters –
perm	– the permutation vector.
nsize	– number of supernodes

◆ MmdUpdate()

static void MmdUpdate	(	int	ehead,
		int	neqns,
		int *	xadj,
		int *	adjncy,
		int	delta,
		int *	mdeg,
		int *	head,
		int *	forward,
		int *	backward,
		int *	qsize,
		int *	list,
		int *	marker,
		int	maxint,
		int *	tag
	)

static

MmdUpdate -— multiple minimum degree update purpose – this routine updates the degrees of nodes after a multiple elimination step. input parameters –.

Parameters

ehead	– int the beginning of the list of eliminated nodes (i.e., newly formed elements).
neqns	– int number of equations.
(xadj,adjncy)	–int adjacency structure.
delta	– int tolerance value for multiple elimination.
maxint	– int maximum machine representable (short) integer. updated parameters –
mdeg	– int new minimum degree after degree update.
(head,forward,backward)	– int degree doubly linked structure.
qsize	– int size of supernode.
list,marker	vector for degree update.
tag	– int tag value.

◆ Permute()

int Permute	(	int	N,
		int *	Perm,
		double *	Vector
	)

This function permute a vector, for a given permutation vector and compute the inverse permutation vector. Of course this is a trivial task if we have a second vector to store the new permuted vector, but we do not want to allocated extra memory so that we use a somewhat expensive algorithm which do not yet require a second storage vector. Here we assume that the dimension of permutation array is less than 2**31, because we use the 31th bit of each index of the permutation vector to store a flag, which tells us that we have already permuted that element. The algorithm is very simple, we take one element, look if we have already permuted that element, and if not we permute all elements connected to this element by a permutation cycle i.e. we shift right by one all elements of that cycle. ATTENTION: This function not only permute the given vector, but also compute the inverse permutation vector which is stored at the place of the permutation vector.

Parameters

N	int
Perm	int
Vector	double

Returns: int

◆ ReleaseBlockMatrix()

int ReleaseBlockMatrix ( SD_BLOCK_MATRIX_DATA * pMat )

Release all the data allocated for the numerical factorization algorithm.

Parameters

pMat	SD_BLOCK_MATRIX_DATA

Returns: int

◆ ReleaseConMatrix()

int ReleaseConMatrix ( SD_CON_MATRIX_DATA * pMat )

◆ RunMmd()

static void RunMmd	(	int	neqns,
		int *	xadj,
		int *	adjncy,
		int *	invp,
		int *	perm,
		int	delta,
		int *	head,
		int *	qsize,
		int *	nsize,
		int *	list,
		int *	marker,
		int	maxint,
		int *	ncsub
	)

static

RunMmd – multiple minimum external degree purpose – this routine implements the minimum degree algorithm. it makes use of the implicit representation of elimination graphs by quotient graphs, and the notion of indistinguishable nodes. It also implements the modifications by multiple elimination and minimum external degree. Caution – the adjacency vector adjncy will be destroyed. Input parameters –.

Parameters

neqns	– number of equations.
(xadj,adjncy)	– the adjacency structure.
delta	– tolerance value for multiple elimination.
maxint	– maximum machine representable (short) integer (any smaller estimate will do) for marking nodes. Output parameters –
nsize	– number of supernodes.
perm	– the minimum degree ordering – used temporarily for degree backward link.
invp	– the inverse of perm – used temporarily for degree forward link.
ncsub	– an upper bound on the number of nonzero subscripts for the compressed storage scheme. Working parameters –
head	– vector for head of degree lists.
qsize	– vector for size of supernodes.
list	– vector for temporary linked lists.
marker	– a temporary marker vector. Subroutines used – MmdElimin, MmdInit, MmdNumbering, MmdUpdate.

◆ SymbolicFact()

int SymbolicFact ( SD_MATRIX_DATA * pMat )

inline

Variable Documentation

◆ gd_MemErr

bool gd_MemErr

static

Classes

Macros

Functions

Variables

Macro Definition Documentation

◆ BLOCK_C0

◆ BLOCK_C1

◆ BLOCK_INIT

◆ BLOCK_JUMP

◆ BLOCK_NEXT

◆ BLOCK_SIZE

◆ DIAGONAL

◆ ERROR_SOLVER

◆ EXIT

◆ FACT_SYM_MAT

◆ FIND_COL_BLOCK

◆ FIRST_BLOCK_ROW

◆ GD_FREE

◆ GD_MALLOC

◆ GD_REALLOC

◆ LAST_BLOCK_ROW

◆ pC0_FIRST_COL

◆ pFIRST_COL_BLOCK

◆ pSIZE_FIRST_COL

◆ SD_ALLOC_CHUNK

◆ SD_ALLOC_COL

◆ SD_ALLOC_COL_BLOCK

◆ SD_CHUNK_REALLOC_SIZE

◆ SD_COL

◆ SD_DESTROY_CHUNK

◆ SD_FIND_COL

◆ SD_FREE_COL_BLOCK

◆ SD_FREE_COL_BLOCK_0

◆ SD_GET_COL

◆ SD_GET_COL_BLOCK

◆ SD_INSERT_COL

◆ SD_MARKED

◆ SD_N_ALLOC_COL

◆ SD_N_ALLOC_COL_BLOCK

◆ SD_P_FIRST_COL_BLOCK

◆ SD_P_SIZE_COL_BLOCK

◆ SD_ROW

◆ SD_SEARCH_BLOCK_ROW

◆ SEARCH_COL

◆ SEARCH_ROW

◆ USER_ERROR

◆ VD_AXPY

◆ VD_AXPY_JUMP

◆ VD_AXPY_POS

◆ VD_DOT_POS

Function Documentation

◆ AllocateConData()

◆ BuildSparseConFormat()

◆ ComputeBlockMatrix()

◆ ComputeFillIn()

◆ ComputePermutation()

◆ ComputeTmpConMatrix()

◆ ds_AssembleMatrix()

◆ ds_DefineConnectivity()

◆ ds_Initialize()

◆ ds_Solve()

◆ FACT_SYM_MAT_BLOCK()

◆ InverseMatrixVector()

◆ InvertMatrix()

◆ MERGE_COL_BLOCK()

◆ MmdElimin()

◆ MmdInit()

◆ MmdNumbering()

◆ MmdUpdate()

◆ Permute()

◆ ReleaseBlockMatrix()

◆ ReleaseConMatrix()

◆ RunMmd()

◆ SymbolicFact()

Variable Documentation

◆ gd_MemErr