SNOWPACK 3.7.0
Classes | Macros | Functions | Variables
Solver.cc File Reference
#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 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. 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; \
} \
} \
}
VD_AXPY
#define VD_AXPY(N, A, X, Y)
Definition: Solver.cc:93

◆ 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_CHUNK_REALLOC_SIZE
#define SD_CHUNK_REALLOC_SIZE
Definition: Solver.cc:121

◆ 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_CHUNK
#define SD_ALLOC_CHUNK(CHUNK, SIZE)
Definition: Solver.cc:123
SD_COL_DATA
Definition: Solver.h:91

◆ 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_FREE_COL_BLOCK_0
#define SD_FREE_COL_BLOCK_0(pCOL_BLOCK, pMAT)
Definition: Solver.cc:208
SD_COL_BLOCK_DATA
Definition: Solver.h:117

◆ 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; \
}
GD_FREE
#define GD_FREE(pointer)
Definition: Solver.cc:74

◆ 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; } \
} \
}
SD_COL
#define SD_COL(pCOL)
Definition: Solver.cc:182
SD_COL_DATA::Next
struct SD_COL_DATA * Next
Definition: Solver.h:93

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_N_ALLOC_COL
#define SD_N_ALLOC_COL
Definition: Solver.cc:187
SD_ALLOC_COL
#define SD_ALLOC_COL(N_COL, pMAT)
Definition: Solver.cc:189

◆ 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_N_ALLOC_COL_BLOCK
#define SD_N_ALLOC_COL_BLOCK
Definition: Solver.cc:222
SD_ALLOC_COL_BLOCK
#define SD_ALLOC_COL_BLOCK(N_COL_BLOCK, pMAT)
Definition: Solver.cc:214

◆ 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; } \
} \
}
SD_ROW_BLOCK_DATA
The data structure to store the matrix for numerical factorization is a simple one....
Definition: Solver.h:66
SD_ROW_BLOCK_DATA::Row1
int Row1
Definition: Solver.h:68

◆ 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_; } \
} \
}
SD_P_FIRST_COL_BLOCK
#define SD_P_FIRST_COL_BLOCK(pMAT, pROW)
Definition: Solver.cc:163
DIAGONAL
#define DIAGONAL(DIM, K)
Definition: Solver.cc:358
SD_P_SIZE_COL_BLOCK
#define SD_P_SIZE_COL_BLOCK(pMAT, pROW)
Definition: Solver.cc:164

◆ 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
Dimint
pMatSD_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
pMatSD_CON_MATRIX_DATA
pRowStart0int
pColumn0int
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
pTmpMatSD_TMP_CON_MATRIX_DATA
pMatSD_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
pMatSD_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
pMatSD_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
pMat0SD_CON_MATRIX_DATA
pMatSD_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]pMat0pointer to the matrix [A] opaque data returned by ds_Initialize()
[in]nEqno. of equations for one element forming a crique
[in]EqElement list of equations for one element.
[in]Dimfirst dimension of the 2D-array ElMat[][Dim]
[in]ElMatelement 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]pMat0pointer to the matrix [A] opaque data returned by ds_Initialize()
[in]nEqNo. of equations for one element forming a crique
[in]EqElement list of equations for more elements with equal no. of eqs.
[in]nElNo. of elements ( 0 <= i "<" nEq ; 0 <= e "<" nEl )
[in]Dimfirst 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
MatDimdimension of the matrix [A]. The true number of equations and unknowns is given by: MatDim * Multiplicity
ppMatA 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]Codefunctionlity code defined above
[in]pMatpointer to the matrix [A] opaque data
[in]pXright 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
pMatSD_BLOCK_MATRIX_DATA
Xdouble
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
pMatSD_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,markervector 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
Nint
Permint
Vectordouble
Returns
int

◆ ReleaseBlockMatrix()

int ReleaseBlockMatrix ( SD_BLOCK_MATRIX_DATA *  pMat)

Release all the data allocated for the numerical factorization algorithm.

Parameters
pMatSD_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