luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
ngspice/src/maths/KLU/klusmp.c

1805 lines
59 KiB
C
Raw Blame History

/* KLU to SMP Interface
* Francesco Lannutti
* July 2020
*/
#include "ngspice/config.h"
#include <assert.h>
#include <stdio.h>
#include <math.h>
#include "ngspice/spmatrix.h"
#include "../sparse/spdefs.h"
#include "ngspice/smpdefs.h"
#if defined (_MSC_VER)
extern double scalbn(double, int);
#define logb _logb
extern double logb(double);
#endif
static void LoadGmin_CSC (double **diag, unsigned int n, double Gmin) ;
static void LoadGmin (SMPmatrix *eMatrix, double Gmin) ;
typedef struct sElement {
unsigned int row ;
unsigned int col ;
double *pointer ;
unsigned int group ;
} Element ;
static int
CompareRow (const void *a, const void *b)
{
Element *A = (Element *) a ;
Element *B = (Element *) b ;
return
(A->row > B->row) ? 1 :
(A->row < B->row) ? -1 :
0 ;
}
static int
CompareColumn (const void *a, const void *b)
{
Element *A = (Element *) a ;
Element *B = (Element *) b ;
return
(A->col > B->col) ? 1 :
(A->col < B->col) ? -1 :
0 ;
}
static void
Compress (unsigned int *Ai, unsigned int *Bp, unsigned int n, unsigned int nz)
{
unsigned int i, j ;
for (i = 0 ; i <= Ai [0] ; i++)
Bp [i] = 0 ;
j = Ai [0] + 1 ;
for (i = 1 ; i < nz ; i++)
{
if (Ai [i] == Ai [i - 1] + 1)
{
Bp [j] = i ;
j++ ;
}
else if (Ai [i] > Ai [i - 1] + 1)
{
for ( ; j <= Ai [i] ; j++)
Bp [j] = i ;
}
}
for ( ; j <= n ; j++)
Bp [j] = i ;
}
int
BindCompare (const void *a, const void *b)
{
BindElement *A = (BindElement *) a ;
BindElement *B = (BindElement *) b ;
return
(A->COO > B->COO) ? 1 :
(A->COO < B->COO) ? -1 :
0 ;
}
#ifdef CIDER
int
BindCompareKLUforCIDER (const void *a, const void *b)
{
BindElementKLUforCIDER *A = (BindElementKLUforCIDER *) a ;
BindElementKLUforCIDER *B = (BindElementKLUforCIDER *) b ;
return
(A->COO > B->COO) ? 1 :
(A->COO < B->COO) ? -1 :
0 ;
}
int
BindKluCompareCSCKLUforCIDER (const void *a, const void *b)
{
BindElementKLUforCIDER *A = (BindElementKLUforCIDER *) a ;
BindElementKLUforCIDER *B = (BindElementKLUforCIDER *) b ;
return
(A->CSC_Complex > B->CSC_Complex) ? 1 :
(A->CSC_Complex < B->CSC_Complex) ? -1 :
0 ;
}
#endif
void SMPconvertCOOtoCSC (SMPmatrix *Matrix)
{
Element *MatrixCOO ;
KluLinkedListCOO *current, *temp ;
unsigned int *Ap_COO, current_group, i, j ;
/* Allocate the compressed COO elements */
MatrixCOO = (Element *) malloc (Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ * sizeof (Element)) ;
/* Populate the compressed COO elements and COO value of Binding Table */
/* Delete the Linked List in the meantime */
i = 0 ;
temp = Matrix->SMPkluMatrix->KLUmatrixLinkedListCOO ;
while (temp != NULL) {
MatrixCOO [i].row = temp->row ;
MatrixCOO [i].col = temp->col ;
MatrixCOO [i].pointer = temp->pointer ;
MatrixCOO [i].group = 0 ;
current = temp ;
temp = temp->next ;
free (current) ;
current = NULL ;
i++ ;
}
/* Order the MatrixCOO along the columns */
qsort (MatrixCOO, Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ, sizeof (Element), CompareColumn) ;
/* Order the MatrixCOO along the rows */
i = 0 ;
while (i < Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ)
{
/* Look for the next column */
for (j = i + 1 ; j < Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ ; j++)
{
if (MatrixCOO [j].col != MatrixCOO [i].col)
{
break ;
}
}
qsort (MatrixCOO + i, j - i, sizeof (Element), CompareRow) ;
i = j ;
}
/* Assign labels to avoid duplicates */
for (i = 0, j = 1 ; i < Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ - 1 ; i++, j++) {
if ((MatrixCOO [i].col == MatrixCOO [j].col) && (MatrixCOO [i].row == MatrixCOO [j].row)) {
// If col and row are the same
MatrixCOO [j].group = MatrixCOO [i].group ;
} else if ((MatrixCOO [i].col != MatrixCOO [j].col) || (MatrixCOO [i].row != MatrixCOO [j].row)) {
// If or col either row are different, it isn't a duplicate, so assign the next label and store it in 'nz'
MatrixCOO [j].group = MatrixCOO [i].group + 1 ;
} else {
printf ("Error: Strange behavior during label assignment\n") ;
}
}
/* Assign N and NZ */
Matrix->SMPkluMatrix->KLUmatrixN = (unsigned int)MatrixCOO [Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ - 1].col + 1 ;
Matrix->SMPkluMatrix->KLUmatrixNZ = MatrixCOO [Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ - 1].group + 1 ;
/* Allocate Diag Gmin CSC Vector */
Matrix->SMPkluMatrix->KLUmatrixDiag = (double **) malloc (Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double *)) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
Matrix->SMPkluMatrix->KLUmatrixDiag [i] = NULL ;
}
/* Allocate the temporary COO Column Index */
Ap_COO = (unsigned int *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof (unsigned int)) ;
/* Allocate the needed KLU data structures */
Matrix->SMPkluMatrix->KLUmatrixAp = (int *) malloc ((Matrix->SMPkluMatrix->KLUmatrixN + 1) * sizeof (int)) ;
Matrix->SMPkluMatrix->KLUmatrixAi = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof (int)) ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = (BindElement *) malloc (Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ * sizeof (BindElement)) ;
Matrix->SMPkluMatrix->KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof (double)) ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof (double)) ;
Matrix->SMPkluMatrix->KLUmatrixIntermediate = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double)) ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double)) ;
/* Copy back the Matrix in partial CSC */
for (i = 0, current_group = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ ; i++)
{
if (MatrixCOO [i].group > current_group) {
current_group = MatrixCOO [i].group ;
}
Ap_COO [current_group] = MatrixCOO [i].col ;
Matrix->SMPkluMatrix->KLUmatrixAi [current_group] = (int)MatrixCOO [i].row ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].COO = MatrixCOO [i].pointer ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].CSC = &(Matrix->SMPkluMatrix->KLUmatrixAx [current_group]) ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].CSC_Complex = &(Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * current_group]) ;
if (MatrixCOO [i].col == MatrixCOO [i].row) {
Matrix->SMPkluMatrix->KLUmatrixDiag [MatrixCOO [i].col] = Matrix->SMPkluMatrix->KLUmatrixBindStructCOO [i].CSC ;
}
}
/* Compress the COO Column Index to CSC Column Index */
Compress (Ap_COO, (unsigned int *)Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixN, Matrix->SMPkluMatrix->KLUmatrixNZ) ;
/* Free the temporary stuff */
free (Ap_COO) ;
free (MatrixCOO) ;
/* Sort the Binding Table */
qsort (Matrix->SMPkluMatrix->KLUmatrixBindStructCOO, Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ, sizeof (BindElement), BindCompare) ;
/* Set the Matrix as Real */
Matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal ;
return ;
}
#ifdef CIDER
typedef struct sElementKLUforCIDER {
unsigned int row ;
unsigned int col ;
double *pointer ;
} ElementKLUforCIDER ;
void SMPconvertCOOtoCSCKLUforCIDER (SMPmatrix *Matrix)
{
ElementKLUforCIDER *MatrixCOO ;
unsigned int *Ap_COO, i, j, nz ;
/* Count the non-zero elements and store it */
nz = 0 ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
if ((Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER [i] != -1) && (Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER [i] != -1)) {
nz++ ;
}
}
Matrix->SMPkluMatrix->KLUmatrixNZ = nz ;
/* Allocate the compressed COO elements */
MatrixCOO = (ElementKLUforCIDER *) malloc (nz * sizeof (ElementKLUforCIDER)) ;
/* Allocate the temporary COO Column Index */
Ap_COO = (unsigned int *) malloc (nz * sizeof (unsigned int)) ;
/* Allocate the needed KLU data structures */
Matrix->SMPkluMatrix->KLUmatrixAp = (int *) malloc ((Matrix->SMPkluMatrix->KLUmatrixN + 1) * sizeof (int)) ;
Matrix->SMPkluMatrix->KLUmatrixAi = (int *) malloc (nz * sizeof (int)) ;
Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER = (BindElementKLUforCIDER *) malloc (nz * sizeof (BindElementKLUforCIDER)) ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = (double *) malloc (2 * nz * sizeof (double)) ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double)) ;
/* Populate the compressed COO elements and COO value of Binding Table */
j = 0 ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
if ((Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER [i] != -1) && (Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER [i] != -1)) {
MatrixCOO [j].row = (unsigned int)Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER [i] ;
MatrixCOO [j].col = (unsigned int)Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER [i] ;
MatrixCOO [j].pointer = &(Matrix->SMPkluMatrix->KLUmatrixValueComplexCOOforCIDER [2 * i]) ;
j++ ;
}
}
/* Order the MatrixCOO along the columns */
qsort (MatrixCOO, nz, sizeof (ElementKLUforCIDER), CompareColumn) ;
/* Order the MatrixCOO along the rows */
i = 0 ;
while (i < nz)
{
/* Look for the next column */
for (j = i + 1 ; j < nz ; j++)
{
if (MatrixCOO [j].col != MatrixCOO [i].col)
{
break ;
}
}
qsort (MatrixCOO + i, j - i, sizeof (ElementKLUforCIDER), CompareRow) ;
i = j ;
}
/* Copy back the Matrix in partial CSC */
for (i = 0 ; i < nz ; i++)
{
Ap_COO [i] = MatrixCOO [i].col ;
Matrix->SMPkluMatrix->KLUmatrixAi [i] = (int)MatrixCOO [i].row ;
Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER [i].COO = MatrixCOO [i].pointer ;
Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER [i].CSC_Complex = &(Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i]) ;
}
/* Compress the COO Column Index to CSC Column Index */
Compress (Ap_COO, (unsigned int *)Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixN, nz) ;
/* Free the temporary stuff */
free (Ap_COO) ;
free (MatrixCOO) ;
/* Sort the Binding Table */
qsort (Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER, nz, sizeof (BindElementKLUforCIDER), BindCompareKLUforCIDER) ;
return ;
}
#endif
/*
* SMPmakeElt()
*/
double *
SMPmakeElt (SMPmatrix *Matrix, int Row, int Col)
{
KluLinkedListCOO *temp ;
if (Matrix->CKTkluMODE) {
if ((Row > 0) && (Col > 0)) {
Row = Row - 1 ;
Col = Col - 1 ;
temp = (KluLinkedListCOO *) malloc (sizeof (KluLinkedListCOO)) ;
temp->row = (unsigned int)Row ;
temp->col = (unsigned int)Col ;
temp->pointer = (double *) malloc (sizeof (double)) ;
temp->next = Matrix->SMPkluMatrix->KLUmatrixLinkedListCOO ;
Matrix->SMPkluMatrix->KLUmatrixLinkedListCOO = temp ;
Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ++ ;
return temp->pointer ;
} else {
return Matrix->SMPkluMatrix->KLUmatrixTrashCOO ;
}
} else {
return spGetElement (Matrix->SPmatrix, Row, Col) ;
}
}
#ifdef CIDER
double *
SMPmakeEltKLUforCIDER (SMPmatrix *Matrix, int Row, int Col)
{
if (Matrix->CKTkluMODE) {
if ((Row > 0) && (Col > 0)) {
Row = Row - 1 ;
Col = Col - 1 ;
Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER [Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col] = Row ;
Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER [Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col] = Col ;
return &(Matrix->SMPkluMatrix->KLUmatrixValueComplexCOOforCIDER [2 * (Row * (int)Matrix->SMPkluMatrix->KLUmatrixN + Col)]) ;
} else {
return Matrix->SMPkluMatrix->KLUmatrixTrashCOO ;
}
} else {
return spGetElement (Matrix->SPmatrix, Row, Col) ;
}
}
#endif
/*
* SMPcClear()
*/
void
SMPcClear (SMPmatrix *Matrix)
{
unsigned int i ;
if (Matrix->CKTkluMODE)
{
for (i = 0 ; i < 2 * Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
Matrix->SMPkluMatrix->KLUmatrixAxComplex [i] = 0 ;
}
} else {
spClear (Matrix->SPmatrix) ;
}
}
/*
* SMPclear()
*/
void
SMPclear (SMPmatrix *Matrix)
{
unsigned int i ;
if (Matrix->CKTkluMODE)
{
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
Matrix->SMPkluMatrix->KLUmatrixAx [i] = 0 ;
}
} else {
spClear (Matrix->SPmatrix) ;
}
}
#ifdef CIDER
void
SMPclearKLUforCIDER (SMPmatrix *Matrix)
{
unsigned int i ;
for (i = 0 ; i < 2 * Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
Matrix->SMPkluMatrix->KLUmatrixAxComplex [i] = 0 ;
}
}
#endif
#define NG_IGNORE(x) (void)x
/*
* SMPcLUfac()
*/
/*ARGSUSED*/
int
SMPcLUfac (SMPmatrix *Matrix, double PivTol)
{
int ret ;
NG_IGNORE (PivTol) ;
if (Matrix->CKTkluMODE)
{
ret = klu_z_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (ret == 0)
{
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (ReFactor): KLU Matrix is SINGULAR\n") ;
return E_SINGULAR ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (ReFactor): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
{
fprintf (stderr, "Error (ReFactor): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) {
fprintf (stderr, "Error (ReFactor): KLUnumeric object is NULL. A problem occurred\n") ;
}
return 1 ;
} else {
return 0 ;
}
} else {
spSetComplex (Matrix->SPmatrix) ;
return spFactor (Matrix->SPmatrix) ;
}
}
/*
* SMPluFac()
*/
/*ARGSUSED*/
int
SMPluFac (SMPmatrix *Matrix, double PivTol, double Gmin)
{
int ret ;
NG_IGNORE (PivTol) ;
if (Matrix->CKTkluMODE)
{
LoadGmin_CSC (Matrix->SMPkluMatrix->KLUmatrixDiag, Matrix->SMPkluMatrix->KLUmatrixN, Gmin) ;
ret = klu_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAx,
Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (ret == 0)
{
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (ReFactor): KLU Matrix is SINGULAR\n") ;
return E_SINGULAR ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (ReFactor): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
{
fprintf (stderr, "Error (ReFactor): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) {
fprintf (stderr, "Error (ReFactor): KLUnumeric object is NULL. A problem occurred\n") ;
}
return 1 ;
} else {
return 0 ;
}
} else {
spSetReal (Matrix->SPmatrix) ;
LoadGmin (Matrix, Gmin) ;
return spFactor (Matrix->SPmatrix) ;
}
}
#ifdef CIDER
int
SMPluFacKLUforCIDER (SMPmatrix *Matrix)
{
unsigned int i ;
double *KLUmatrixAx ;
int ret ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
ret = klu_z_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
} else {
/* Allocate the Real Matrix */
KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
/* Copy the Complex Matrix into the Real Matrix */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
}
/* Re-Factor the Real Matrix */
ret = klu_refactor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, KLUmatrixAx,
Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/* Free the Real Matrix Storage */
free (KLUmatrixAx) ;
}
if (ret == 0)
{
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (ReFactor): KLU Matrix is SINGULAR\n") ;
return E_SINGULAR ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (ReFactor): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
{
fprintf (stderr, "Error (ReFactor): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) {
fprintf (stderr, "Error (ReFactor): KLUnumeric object is NULL. A problem occurred\n") ;
}
return 1 ;
} else {
return 0 ;
}
} else {
return spFactor (Matrix->SPmatrix) ;
}
}
#endif
/*
* SMPcReorder()
*/
int
SMPcReorder (SMPmatrix *Matrix, double PivTol, double PivRel, int *NumSwaps)
{
if (Matrix->CKTkluMODE)
{
Matrix->SMPkluMatrix->KLUmatrixCommon->tol = PivTol ;
if (Matrix->SMPkluMatrix->KLUmatrixNumeric != NULL) {
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
}
Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_z_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
Matrix->SMPkluMatrix->KLUmatrixAxComplex, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
Matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL)
{
fprintf (stderr, "Error (Factor): KLUnumeric object is NULL. A problem occurred\n") ;
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (Factor): KLU Matrix is SINGULAR\n") ;
return 0 ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (Factor): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
fprintf (stderr, "Error (Factor): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
fprintf (stderr, "Error (Factor): KLUsymbolic object is NULL. A problem occurred\n") ;
}
return 1 ;
} else {
return 0 ;
}
} else {
*NumSwaps = 1 ;
spSetComplex (Matrix->SPmatrix) ;
return spOrderAndFactor (Matrix->SPmatrix, NULL, (spREAL)PivRel, (spREAL)PivTol, YES) ;
}
}
/*
* SMPreorder()
*/
int
SMPreorder (SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin)
{
if (Matrix->CKTkluMODE)
{
LoadGmin_CSC (Matrix->SMPkluMatrix->KLUmatrixDiag, Matrix->SMPkluMatrix->KLUmatrixN, Gmin) ;
Matrix->SMPkluMatrix->KLUmatrixCommon->tol = PivTol ;
if (Matrix->SMPkluMatrix->KLUmatrixNumeric != NULL) {
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
}
Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
Matrix->SMPkluMatrix->KLUmatrixAx, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
Matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL)
{
fprintf (stderr, "Error (Factor): KLUnumeric object is NULL. A problem occurred\n") ;
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (Factor): KLU Matrix is SINGULAR\n") ;
return 0 ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (Factor): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
fprintf (stderr, "Error (Factor): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
fprintf (stderr, "Error (Factor): KLUsymbolic object is NULL. A problem occurred\n") ;
}
return 1 ;
} else {
return 0 ;
}
} else {
spSetReal (Matrix->SPmatrix) ;
LoadGmin (Matrix, Gmin) ;
return spOrderAndFactor (Matrix->SPmatrix, NULL, (spREAL)PivRel, (spREAL)PivTol, YES) ;
}
}
#ifdef CIDER
int
SMPreorderKLUforCIDER (SMPmatrix *Matrix)
{
unsigned int i ;
double *KLUmatrixAx ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixNumeric != NULL) {
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
}
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_z_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
Matrix->SMPkluMatrix->KLUmatrixAxComplex, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
Matrix->SMPkluMatrix->KLUmatrixCommon) ;
} else {
/* Allocate the Real Matrix */
KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
/* Copy the Complex Matrix into the Real Matrix */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
}
/* Factor the Real Matrix */
Matrix->SMPkluMatrix->KLUmatrixNumeric = klu_factor (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi,
KLUmatrixAx, Matrix->SMPkluMatrix->KLUmatrixSymbolic,
Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/* Free the Real Matrix Storage */
free (KLUmatrixAx) ;
}
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL)
{
fprintf (stderr, "Error (Factor): KLUnumeric object is NULL. A problem occurred\n") ;
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (Factor): KLU Matrix is SINGULAR\n") ;
return 0 ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (Factor): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX) {
fprintf (stderr, "Error (Factor): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
fprintf (stderr, "Error (Factor): KLUsymbolic object is NULL. A problem occurred\n") ;
}
return 1 ;
} else {
return 0 ;
}
} else {
return spFactor (Matrix->SPmatrix) ;
}
}
#endif
/*
* SMPcaSolve()
*/
void
SMPcaSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], double iSpare[])
{
int ret ;
unsigned int i ;
NG_IGNORE (iSpare) ;
NG_IGNORE (Spare) ;
if (Matrix->CKTkluMODE)
{
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] = RHS [i + 1] ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] = iRHS [i + 1] ;
}
ret = klu_z_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
RHS [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] ;
iRHS [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] ;
}
} else {
spSolveTransposed (Matrix->SPmatrix, RHS, RHS, iRHS, iRHS) ;
}
}
/*
* SMPcSolve()
*/
void
SMPcSolve (SMPmatrix *Matrix, double RHS[], double iRHS[], double Spare[], double iSpare[])
{
int ret ;
unsigned int i ;
NG_IGNORE (iSpare) ;
NG_IGNORE (Spare) ;
if (Matrix->CKTkluMODE)
{
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] = RHS [i + 1] ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] = iRHS [i + 1] ;
}
ret = klu_z_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
RHS [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] ;
iRHS [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] ;
}
} else {
spSolve (Matrix->SPmatrix, RHS, RHS, iRHS, iRHS) ;
}
}
/*
* SMPsolve()
*/
void
SMPsolve (SMPmatrix *Matrix, double RHS[], double Spare[])
{
int ret ;
unsigned int i ;
NG_IGNORE (Spare) ;
if (Matrix->CKTkluMODE)
{
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
Matrix->SMPkluMatrix->KLUmatrixIntermediate [i] = RHS [i + 1] ;
}
ret = klu_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
Matrix->SMPkluMatrix->KLUmatrixIntermediate, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (ret == 0)
{
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR) {
fprintf (stderr, "Warning (Solve): KLU Matrix is SINGULAR\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon == NULL) {
fprintf (stderr, "Error (Solve): KLUcommon object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
{
fprintf (stderr, "Error (Solve): KLU Matrix is empty\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixNumeric == NULL) {
fprintf (stderr, "Error (Solve): KLUnumeric object is NULL. A problem occurred\n") ;
}
if (Matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL) {
fprintf (stderr, "Error (Solve): KLUsymbolic object is NULL. A problem occurred\n") ;
}
}
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
RHS [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediate [i] ;
}
} else {
spSolve (Matrix->SPmatrix, RHS, RHS, NULL, NULL) ;
}
}
#ifdef CIDER
void
SMPsolveKLUforCIDER (SMPmatrix *Matrix, double RHS[], double RHSsolution[], double iRHS[], double iRHSsolution[])
{
int ret ;
unsigned int i ;
double *KLUmatrixIntermediate ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex) {
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] = RHS [i + 1] ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] = iRHS [i + 1] ;
}
ret = klu_z_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++)
{
RHSsolution [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i] ;
iRHSsolution [i + 1] = Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex [2 * i + 1] ;
}
} else {
/* Allocate the Intermediate Vector */
KLUmatrixIntermediate = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
KLUmatrixIntermediate [i] = RHS [i + 1] ;
}
ret = klu_solve (Matrix->SMPkluMatrix->KLUmatrixSymbolic, Matrix->SMPkluMatrix->KLUmatrixNumeric, (int)Matrix->SMPkluMatrix->KLUmatrixN, 1,
KLUmatrixIntermediate, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
RHSsolution [i + 1] = KLUmatrixIntermediate [i] ;
}
/* Free the Intermediate Vector */
free (KLUmatrixIntermediate) ;
}
} else {
spSolve (Matrix->SPmatrix, RHS, RHSsolution, iRHS, iRHSsolution) ;
}
}
#endif
/*
* SMPmatSize()
*/
int
SMPmatSize (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE) {
return (int)Matrix->SMPkluMatrix->KLUmatrixN ;
} else {
return spGetSize (Matrix->SPmatrix, 1) ;
}
}
/*
* SMPnewMatrix()
*/
int
SMPnewMatrix (SMPmatrix *Matrix, int size)
{
int Error ;
if (Matrix->CKTkluMODE) {
/* Allocate the KLU Matrix Data Structure */
Matrix->SMPkluMatrix = (KLUmatrix *) malloc (sizeof (KLUmatrix)) ;
Matrix->SMPkluMatrix->KLUmatrixLinkedListNZ = 0 ;
Matrix->SMPkluMatrix->KLUmatrixLinkedListCOO = NULL ;
/* Initialize the KLU Matrix Internal Pointers */
Matrix->SMPkluMatrix->KLUmatrixCommon = (klu_common *) malloc (sizeof (klu_common)) ; ;
Matrix->SMPkluMatrix->KLUmatrixSymbolic = NULL ;
Matrix->SMPkluMatrix->KLUmatrixNumeric = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAx = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal ;
Matrix->SMPkluMatrix->KLUmatrixIntermediate = NULL ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixNZ = 0 ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixDiag = NULL ;
/* Initialize the KLU Common Data Structure */
klu_defaults (Matrix->SMPkluMatrix->KLUmatrixCommon) ;
Matrix->SMPkluMatrix->KLUmatrixCommon->memgrow = Matrix->CKTkluMemGrowFactor ;
/* Allocate KLU data structures */
Matrix->SMPkluMatrix->KLUmatrixN = (unsigned int)size ;
Matrix->SMPkluMatrix->KLUmatrixTrashCOO = (double *) malloc (sizeof (double)) ;
return spOKAY ;
} else {
Matrix->SPmatrix = spCreate (size, 1, &Error) ;
return Error ;
}
}
#ifdef CIDER
int
SMPnewMatrixKLUforCIDER (SMPmatrix *Matrix, int size, unsigned int KLUmatrixIsComplex)
{
int Error ;
unsigned int i ;
if (Matrix->CKTkluMODE) {
/* Allocate the KLU Matrix Data Structure */
Matrix->SMPkluMatrix = (KLUmatrix *) malloc (sizeof (KLUmatrix)) ;
/* Initialize the KLU Matrix Internal Pointers */
Matrix->SMPkluMatrix->KLUmatrixCommon = (klu_common *) malloc (sizeof (klu_common)) ; ;
Matrix->SMPkluMatrix->KLUmatrixSymbolic = NULL ;
Matrix->SMPkluMatrix->KLUmatrixNumeric = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
if (KLUmatrixIsComplex) {
Matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex ;
} else {
Matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUmatrixReal ;
}
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixNZ = 0 ;
Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER = NULL ;
Matrix->SMPkluMatrix->KLUmatrixValueComplexCOOforCIDER = NULL ;
/* Initialize the KLU Common Data Structure */
klu_defaults (Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/* Allocate KLU data structures */
Matrix->SMPkluMatrix->KLUmatrixN = (unsigned int)size ;
Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(int)) ;
Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER = (int *) malloc (Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(int)) ;
Matrix->SMPkluMatrix->KLUmatrixTrashCOO = (double *) malloc (sizeof(double)) ;
Matrix->SMPkluMatrix->KLUmatrixValueComplexCOOforCIDER = (double *) malloc (2 * Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN * sizeof(double)) ;
/* Pre-set the values of Row and Col */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixN * Matrix->SMPkluMatrix->KLUmatrixN ; i++) {
Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER [i] = -1 ;
Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER [i] = -1 ;
}
return spOKAY ;
} else {
Matrix->SPmatrix = spCreate (size, (int)KLUmatrixIsComplex, &Error) ;
return Error ;
}
}
#endif
/*
* SMPdestroy()
*/
void
SMPdestroy (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE)
{
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
klu_free_symbolic (&(Matrix->SMPkluMatrix->KLUmatrixSymbolic), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
free (Matrix->SMPkluMatrix->KLUmatrixAp) ;
free (Matrix->SMPkluMatrix->KLUmatrixAi) ;
free (Matrix->SMPkluMatrix->KLUmatrixAx) ;
free (Matrix->SMPkluMatrix->KLUmatrixAxComplex) ;
free (Matrix->SMPkluMatrix->KLUmatrixIntermediate) ;
free (Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex) ;
free (Matrix->SMPkluMatrix->KLUmatrixBindStructCOO) ;
free (Matrix->SMPkluMatrix->KLUmatrixTrashCOO) ;
Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAx = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixIntermediate = NULL ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixBindStructCOO = NULL ;
Matrix->SMPkluMatrix->KLUmatrixTrashCOO = NULL ;
} else {
spDestroy (Matrix->SPmatrix) ;
}
}
#ifdef CIDER
void
SMPdestroyKLUforCIDER (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE)
{
klu_free_numeric (&(Matrix->SMPkluMatrix->KLUmatrixNumeric), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
klu_free_symbolic (&(Matrix->SMPkluMatrix->KLUmatrixSymbolic), Matrix->SMPkluMatrix->KLUmatrixCommon) ;
free (Matrix->SMPkluMatrix->KLUmatrixAp) ;
free (Matrix->SMPkluMatrix->KLUmatrixAi) ;
free (Matrix->SMPkluMatrix->KLUmatrixAxComplex) ;
free (Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex) ;
free (Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER) ;
free (Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER) ;
free (Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER) ;
free (Matrix->SMPkluMatrix->KLUmatrixValueComplexCOOforCIDER) ;
free (Matrix->SMPkluMatrix->KLUmatrixTrashCOO) ;
Matrix->SMPkluMatrix->KLUmatrixAp = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAi = NULL ;
Matrix->SMPkluMatrix->KLUmatrixAxComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixIntermediateComplex = NULL ;
Matrix->SMPkluMatrix->KLUmatrixBindStructForCIDER = NULL ;
Matrix->SMPkluMatrix->KLUmatrixColCOOforCIDER = NULL ;
Matrix->SMPkluMatrix->KLUmatrixRowCOOforCIDER = NULL ;
Matrix->SMPkluMatrix->KLUmatrixValueComplexCOOforCIDER = NULL ;
Matrix->SMPkluMatrix->KLUmatrixTrashCOO = NULL ;
} else {
spDestroy (Matrix->SPmatrix) ;
}
}
#endif
/*
* SMPpreOrder()
*/
int
SMPpreOrder (SMPmatrix *Matrix)
{
if (Matrix->CKTkluMODE)
{
Matrix->SMPkluMatrix->KLUmatrixSymbolic = klu_analyze ((int)Matrix->SMPkluMatrix->KLUmatrixN, Matrix->SMPkluMatrix->KLUmatrixAp,
Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
if (Matrix->SMPkluMatrix->KLUmatrixSymbolic == NULL)
{
fprintf (stderr, "Error (PreOrder): KLUsymbolic object is NULL. A problem occurred\n") ;
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_EMPTY_MATRIX)
{
return 0 ;
}
return 1 ;
} else {
return 0 ;
}
} else {
spMNA_Preorder (Matrix->SPmatrix) ;
return spError (Matrix->SPmatrix) ;
}
}
/*
* SMPprintRHS()
*/
void
SMPprintRHS (SMPmatrix *Matrix, char *Filename, RealVector RHS, RealVector iRHS)
{
if (!Matrix->CKTkluMODE)
spFileVector (Matrix->SPmatrix, Filename, RHS, iRHS) ;
}
/*
* SMPprint()
*/
void
SMPprint (SMPmatrix *Matrix, char *Filename)
{
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex)
{
klu_z_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
(int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
} else {
klu_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAx,
(int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
}
} else {
if (Filename)
spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ;
else
spPrint (Matrix->SPmatrix, 0, 1, 1) ;
}
}
#ifdef CIDER
void
SMPprintKLUforCIDER (SMPmatrix *Matrix, char *Filename)
{
unsigned int i ;
double *KLUmatrixAx ;
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex)
{
klu_z_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
(int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
} else {
/* Allocate the Real Matrix */
KLUmatrixAx = (double *) malloc (Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof(double)) ;
/* Copy the Complex Matrix into the Real Matrix */
for (i = 0 ; i < Matrix->SMPkluMatrix->KLUmatrixNZ ; i++) {
KLUmatrixAx [i] = Matrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] ;
}
/* Print the Real Matrix */
klu_print (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, KLUmatrixAx, (int)Matrix->SMPkluMatrix->KLUmatrixN, NULL, NULL) ;
/* Free the Real Matrix Storage */
free (KLUmatrixAx) ;
}
} else {
if (Filename)
spFileMatrix (Matrix->SPmatrix, Filename, "Circuit Matrix", 0, 1, 1) ;
else
spPrint (Matrix->SPmatrix, 0, 1, 1) ;
}
}
#endif
/*
* SMPgetError()
*/
void
SMPgetError (SMPmatrix *Matrix, int *Col, int *Row)
{
if (Matrix->CKTkluMODE)
{
*Row = Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1 ;
*Col = Matrix->SMPkluMatrix->KLUmatrixCommon->singular_col + 1 ;
} else {
spWhereSingular (Matrix->SPmatrix, Row, Col) ;
}
}
void
spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant, RealNumber *piDeterminant)
{
int I, Size ;
RealNumber Norm, nr, ni ;
ComplexNumber Pivot, cDeterminant, Udiag ;
int *P, *Q ;
double *Rs, *Ux, *Uz ;
unsigned int nSwap, nSwapP, nSwapQ ;
#define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni))
*pExponent = 0 ;
if (Matrix->SMPkluMatrix->KLUmatrixCommon->status == KLU_SINGULAR)
{
*pDeterminant = 0.0 ;
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex == KLUMatrixComplex)
{
*piDeterminant = 0.0 ;
}
return ;
}
Size = (int)Matrix->SMPkluMatrix->KLUmatrixN ;
I = 0 ;
P = (int *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixN * sizeof (int)) ;
Q = (int *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixN * sizeof (int)) ;
Ux = (double *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double)) ;
Rs = (double *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double)) ;
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex == KLUMatrixComplex) /* Complex Case. */
{
cDeterminant.Real = 1.0 ;
cDeterminant.Imag = 0.0 ;
Uz = (double *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixN * sizeof (double)) ;
/*
int *Lp, *Li, *Up, *Ui, *Fp, *Fi, *P, *Q ;
double *Lx, *Lz, *Ux, *Uz, *Fx, *Fz, *Rs ;
Lp = (int *) malloc (((size_t)Matrix->CKTkluN + 1) * sizeof (int)) ;
Li = (int *) malloc ((size_t)Matrix->CKTkluNumeric->lnz * sizeof (int)) ;
Lx = (double *) malloc ((size_t)Matrix->CKTkluNumeric->lnz * sizeof (double)) ;
Lz = (double *) malloc ((size_t)Matrix->CKTkluNumeric->lnz * sizeof (double)) ;
Up = (int *) malloc (((size_t)Matrix->CKTkluN + 1) * sizeof (int)) ;
Ui = (int *) malloc ((size_t)Matrix->CKTkluNumeric->unz * sizeof (int)) ;
Ux = (double *) malloc ((size_t)Matrix->CKTkluNumeric->unz * sizeof (double)) ;
Uz = (double *) malloc ((size_t)Matrix->CKTkluNumeric->unz * sizeof (double)) ;
Fp = (int *) malloc (((size_t)Matrix->CKTkluN + 1) * sizeof (int)) ;
Fi = (int *) malloc ((size_t)Matrix->CKTkluNumeric->Offp [Matrix->CKTkluN] * sizeof (int)) ;
Fx = (double *) malloc ((size_t)Matrix->CKTkluNumeric->Offp [Matrix->CKTkluN] * sizeof (double)) ;
Fz = (double *) malloc ((size_t)Matrix->CKTkluNumeric->Offp [Matrix->CKTkluN] * sizeof (double)) ;
klu_z_extract (Matrix->CKTkluNumeric, Matrix->CKTkluSymbolic,
Lp, Li, Lx, Lz,
Up, Ui, Ux, Uz,
Fp, Fi, Fx, Fz,
P, Q, Rs, NULL,
Matrix->CKTkluCommon) ;
*/
klu_z_extract_Udiag (Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixSymbolic, Ux, Uz, P, Q, Rs, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
/*
for (I = 0 ; I < Matrix->CKTkluNumeric->lnz ; I++)
{
printf ("L - Value: %-.9g\t%-.9g\n", Lx [I], Lz [I]) ;
}
for (I = 0 ; I < Matrix->CKTkluNumeric->unz ; I++)
{
printf ("U - Value: %-.9g\t%-.9g\n", Ux [I], Uz [I]) ;
}
for (I = 0 ; I < Matrix->CKTkluNumeric->Offp [Matrix->CKTkluN] ; I++)
{
printf ("F - Value: %-.9g\t%-.9g\n", Fx [I], Fz [I]) ;
}
for (I = 0 ; I < Matrix->CKTkluN ; I++)
{
printf ("U - Value: %-.9g\t%-.9g\n", Ux [I], Uz [I]) ;
}
*/
nSwapP = 0 ;
for (I = 0 ; I < (int)Matrix->SMPkluMatrix->KLUmatrixN ; I++)
{
if (P [I] != I)
{
nSwapP++ ;
}
}
nSwapP /= 2 ;
nSwapQ = 0 ;
for (I = 0 ; I < (int)Matrix->SMPkluMatrix->KLUmatrixN ; I++)
{
if (Q [I] != I)
{
nSwapQ++ ;
}
}
nSwapQ /= 2 ;
nSwap = nSwapP + nSwapQ ;
/*
free (Lp) ;
free (Li) ;
free (Lx) ;
free (Lz) ;
free (Up) ;
free (Ui) ;
free (Fp) ;
free (Fi) ;
free (Fx) ;
free (Fz) ;
*/
I = 0 ;
while (I < Size)
{
Udiag.Real = 1 / (Ux [I] * Rs [I]) ;
Udiag.Imag = Uz [I] * Rs [I] ;
// printf ("Udiag.Real: %-.9g\tUdiag.Imag %-.9g\n", Udiag.Real, Udiag.Imag) ;
CMPLX_RECIPROCAL (Pivot, Udiag) ;
CMPLX_MULT_ASSIGN (cDeterminant, Pivot) ;
// printf ("cDeterminant.Real: %-.9g\tcDeterminant.Imag %-.9g\n", cDeterminant.Real, cDeterminant.Imag) ;
/* Scale Determinant. */
Norm = NORM (cDeterminant) ;
if (Norm != 0.0)
{
while (Norm >= 1.0e12)
{
cDeterminant.Real *= 1.0e-12 ;
cDeterminant.Imag *= 1.0e-12 ;
*pExponent += 12 ;
Norm = NORM (cDeterminant) ;
}
while (Norm < 1.0e-12)
{
cDeterminant.Real *= 1.0e12 ;
cDeterminant.Imag *= 1.0e12 ;
*pExponent -= 12 ;
Norm = NORM (cDeterminant) ;
}
}
I++ ;
}
/* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */
Norm = NORM (cDeterminant) ;
if (Norm != 0.0)
{
while (Norm >= 10.0)
{
cDeterminant.Real *= 0.1 ;
cDeterminant.Imag *= 0.1 ;
(*pExponent)++ ;
Norm = NORM (cDeterminant) ;
}
while (Norm < 1.0)
{
cDeterminant.Real *= 10.0 ;
cDeterminant.Imag *= 10.0 ;
(*pExponent)-- ;
Norm = NORM (cDeterminant) ;
}
}
if (nSwap % 2 != 0)
{
CMPLX_NEGATE (cDeterminant) ;
}
*pDeterminant = cDeterminant.Real ;
*piDeterminant = cDeterminant.Imag ;
free (Uz) ;
}
else
{
/* Real Case. */
*pDeterminant = 1.0 ;
klu_extract_Udiag (Matrix->SMPkluMatrix->KLUmatrixNumeric, Matrix->SMPkluMatrix->KLUmatrixSymbolic, Ux, P, Q, Rs, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
nSwapP = 0 ;
for (I = 0 ; I < (int)Matrix->SMPkluMatrix->KLUmatrixN ; I++)
{
if (P [I] != I)
{
nSwapP++ ;
}
}
nSwapP /= 2 ;
nSwapQ = 0 ;
for (I = 0 ; I < (int)Matrix->SMPkluMatrix->KLUmatrixN ; I++)
{
if (Q [I] != I)
{
nSwapQ++ ;
}
}
nSwapQ /= 2 ;
nSwap = nSwapP + nSwapQ ;
while (I < Size)
{
*pDeterminant /= (Ux [I] * Rs [I]) ;
/* Scale Determinant. */
if (*pDeterminant != 0.0)
{
while (ABS(*pDeterminant) >= 1.0e12)
{
*pDeterminant *= 1.0e-12 ;
*pExponent += 12 ;
}
while (ABS(*pDeterminant) < 1.0e-12)
{
*pDeterminant *= 1.0e12 ;
*pExponent -= 12 ;
}
}
I++ ;
}
/* Scale Determinant again, this time to be between 1.0 <= x <
10.0. */
if (*pDeterminant != 0.0)
{
while (ABS(*pDeterminant) >= 10.0)
{
*pDeterminant *= 0.1 ;
(*pExponent)++ ;
}
while (ABS(*pDeterminant) < 1.0)
{
*pDeterminant *= 10.0 ;
(*pExponent)-- ;
}
}
if (nSwap % 2 != 0)
{
*pDeterminant = -*pDeterminant ;
}
}
free (P) ;
free (Q) ;
free (Ux) ;
free (Rs) ;
}
/*
* SMPcProdDiag()
* note: obsolete for Spice3d2 and later
*/
int
SMPcProdDiag (SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
if (Matrix->CKTkluMODE)
{
spDeterminant_KLU (Matrix, pExponent, &(pMantissa->real), &(pMantissa->imag)) ;
} else {
spDeterminant (Matrix->SPmatrix, pExponent, &(pMantissa->real), &(pMantissa->imag)) ;
}
return spError (Matrix->SPmatrix) ;
}
/*
* SMPcDProd()
*/
int
SMPcDProd (SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent)
{
double re, im, x, y, z;
int p;
if (Matrix->CKTkluMODE)
{
spDeterminant_KLU (Matrix, &p, &re, &im) ;
} else {
spDeterminant (Matrix->SPmatrix, &p, &re, &im) ;
}
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#endif
#ifndef M_LN10
#define M_LN10 2.30258509299404568402
#endif
#ifdef debug_print
printf ("Determinant 10: (%20g,%20g)^%d\n", re, im, p) ;
#endif
/* Convert base 10 numbers to base 2 numbers, for comparison */
y = p * M_LN10 / M_LN2;
x = (int) y;
y -= x;
/* ASSERT
* x = integral part of exponent, y = fraction part of exponent
*/
/* Fold in the fractional part */
#ifdef debug_print
printf (" ** base10 -> base2 int = %g, frac = %20g\n", x, y) ;
#endif
z = pow (2.0, y) ;
re *= z ;
im *= z ;
#ifdef debug_print
printf (" ** multiplier = %20g\n", z) ;
#endif
/* Re-normalize (re or im may be > 2.0 or both < 1.0 */
if (re != 0.0)
{
y = logb (re) ;
if (im != 0.0)
z = logb (im) ;
else
z = 0 ;
} else if (im != 0.0) {
z = logb (im) ;
y = 0 ;
} else {
/* Singular */
/*printf("10 -> singular\n");*/
y = 0 ;
z = 0 ;
}
#ifdef debug_print
printf (" ** renormalize changes = %g,%g\n", y, z) ;
#endif
if (y < z)
y = z ;
*pExponent = (int)(x + y) ;
x = scalbn (re, (int) -y) ;
z = scalbn (im, (int) -y) ;
#ifdef debug_print
printf (" ** values are: re %g, im %g, y %g, re' %g, im' %g\n", re, im, y, x, z) ;
#endif
pMantissa->real = scalbn (re, (int) -y) ;
pMantissa->imag = scalbn (im, (int) -y) ;
#ifdef debug_print
printf ("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real, pMantissa->imag, *pExponent) ;
#endif
if (Matrix->CKTkluMODE)
{
return 0 ;
} else {
return spError (Matrix->SPmatrix) ;
}
}
/*
* The following routines need internal knowledge of the Sparse data
* structures.
*/
/*
* LOAD GMIN
*
* This routine adds Gmin to each diagonal element. Because Gmin is
* added to the current diagonal, which may bear little relation to
* what the outside world thinks is a diagonal, and because the
* elements that are diagonals may change after calling spOrderAndFactor,
* use of this routine is not recommended. It is included here simply
* for compatibility with Spice3.
*/
static void
LoadGmin_CSC (double **diag, unsigned int n, double Gmin)
{
unsigned int i ;
if (Gmin != 0.0) {
for (i = 0 ; i < n ; i++) {
if (diag [i] != NULL) {
// Not all the elements on the diagonal are present, when the circuit is parsed
*(diag [i]) += Gmin ;
}
}
}
}
static void
LoadGmin (SMPmatrix *eMatrix, double Gmin)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
int I ;
ArrayOfElementPtrs Diag ;
ElementPtr diag ;
/* Begin `LoadGmin'. */
assert (IS_SPARSE (Matrix)) ;
if (Gmin != 0.0) {
Diag = Matrix->Diag ;
for (I = Matrix->Size ; I > 0 ; I--)
{
if ((diag = Diag [I]) != NULL)
diag->Real += Gmin ;
}
}
return ;
}
/*
* FIND ELEMENT
*
* This routine finds an element in the matrix by row and column number.
* If the element exists, a pointer to it is returned. If not, then NULL
* is returned unless the CreateIfMissing flag is TRUE, in which case a
* pointer to the new element is returned.
*/
SMPelement *
SMPfindElt (SMPmatrix *eMatrix, int Row, int Col, int CreateIfMissing)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
if (eMatrix->CKTkluMODE)
{
int i ;
Row = Row - 1 ;
Col = Col - 1 ;
if ((Row < 0) || (Col < 0)) {
printf ("Error: Cannot find an element with row '%d' and column '%d' in the KLU matrix\n", Row, Col) ;
return NULL ;
}
for (i = eMatrix->SMPkluMatrix->KLUmatrixAp [Col] ; i < eMatrix->SMPkluMatrix->KLUmatrixAp [Col + 1] ; i++) {
if (eMatrix->SMPkluMatrix->KLUmatrixAi [i] == Row) {
if (eMatrix->SMPkluMatrix->KLUmatrixIsComplex == KLUmatrixReal) {
return (SMPelement *) &(eMatrix->SMPkluMatrix->KLUmatrixAx [i]) ;
} else if (eMatrix->SMPkluMatrix->KLUmatrixIsComplex == KLUMatrixComplex) {
return (SMPelement *) &(eMatrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i]) ;
} else {
return NULL ;
}
}
}
return NULL ;
} else {
ElementPtr Element ;
/* Begin `SMPfindElt'. */
assert (IS_SPARSE (Matrix)) ;
Row = Matrix->ExtToIntRowMap [Row] ;
Col = Matrix->ExtToIntColMap [Col] ;
Element = Matrix->FirstInCol [Col] ;
Element = spcFindElementInCol (Matrix, &Element, Row, Col, CreateIfMissing) ;
return (SMPelement *)Element ;
}
}
/* XXX The following should probably be implemented in spUtils */
/*
* SMPcZeroCol()
*/
int
SMPcZeroCol (SMPmatrix *eMatrix, int Col)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Element ;
if (eMatrix->CKTkluMODE)
{
int i ;
for (i = eMatrix->SMPkluMatrix->KLUmatrixAp [Col - 1] ; i < eMatrix->SMPkluMatrix->KLUmatrixAp [Col] ; i++)
{
eMatrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i] = 0 ;
eMatrix->SMPkluMatrix->KLUmatrixAxComplex [2 * i + 1] = 0 ;
}
return 0 ;
} else {
Col = Matrix->ExtToIntColMap [Col] ;
for (Element = Matrix->FirstInCol [Col] ; Element != NULL ; Element = Element->NextInCol)
{
Element->Real = 0.0 ;
Element->Imag = 0.0 ;
}
return spError (Matrix) ;
}
}
/*
* SMPcAddCol()
*/
int
SMPcAddCol (SMPmatrix *eMatrix, int Accum_Col, int Addend_Col)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Accum, Addend, *Prev ;
Accum_Col = Matrix->ExtToIntColMap [Accum_Col] ;
Addend_Col = Matrix->ExtToIntColMap [Addend_Col] ;
Addend = Matrix->FirstInCol [Addend_Col] ;
Prev = &Matrix->FirstInCol [Accum_Col] ;
Accum = *Prev;
while (Addend != NULL)
{
while (Accum && Accum->Row < Addend->Row)
{
Prev = &Accum->NextInCol ;
Accum = *Prev ;
}
if (!Accum || Accum->Row > Addend->Row)
{
Accum = spcCreateElement (Matrix, Addend->Row, Accum_Col, Prev, 0) ;
}
Accum->Real += Addend->Real ;
Accum->Imag += Addend->Imag ;
Addend = Addend->NextInCol ;
}
return spError (Matrix) ;
}
/*
* SMPzeroRow()
*/
int
SMPzeroRow (SMPmatrix *eMatrix, int Row)
{
MatrixPtr Matrix = eMatrix->SPmatrix ;
ElementPtr Element ;
Row = Matrix->ExtToIntColMap [Row] ;
if (Matrix->RowsLinked == NO)
spcLinkRows (Matrix) ;
if (Matrix->PreviousMatrixWasComplex || Matrix->Complex)
{
for (Element = Matrix->FirstInRow[Row] ; Element != NULL; Element = Element->NextInRow)
{
Element->Real = 0.0 ;
Element->Imag = 0.0 ;
}
} else {
for (Element = Matrix->FirstInRow [Row] ; Element != NULL ; Element = Element->NextInRow)
{
Element->Real = 0.0 ;
}
}
return spError (Matrix) ;
}
/*
* SMPconstMult()
*/
void
SMPconstMult (SMPmatrix *Matrix, double constant)
{
if (Matrix->CKTkluMODE)
{
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex)
{
klu_z_constant_multiply (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAxComplex,
(int)Matrix->SMPkluMatrix->KLUmatrixN, Matrix->SMPkluMatrix->KLUmatrixCommon, constant) ;
} else {
klu_constant_multiply (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAx,
(int)Matrix->SMPkluMatrix->KLUmatrixN, Matrix->SMPkluMatrix->KLUmatrixCommon, constant) ;
}
} else {
spConstMult (Matrix->SPmatrix, constant) ;
}
}
/*
* SMPmultiply()
*/
void
SMPmultiply (SMPmatrix *Matrix, double *RHS, double *Solution, double *iRHS, double *iSolution)
{
if (Matrix->CKTkluMODE)
{
int *Ap_CSR, *Ai_CSR ;
double *Ax_CSR ;
Ap_CSR = (int *) malloc ((size_t)(Matrix->SMPkluMatrix->KLUmatrixN + 1) * sizeof (int)) ;
Ai_CSR = (int *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof (int)) ;
if (Matrix->SMPkluMatrix->KLUmatrixIsComplex)
{
Ax_CSR = (double *) malloc ((size_t)(2 * Matrix->SMPkluMatrix->KLUmatrixNZ) * sizeof (double)) ;
klu_z_convert_matrix_in_CSR (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAxComplex, Ap_CSR,
Ai_CSR, Ax_CSR, (int)Matrix->SMPkluMatrix->KLUmatrixN, (int)Matrix->SMPkluMatrix->KLUmatrixNZ, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
klu_z_matrix_vector_multiply (Ap_CSR, Ai_CSR, Ax_CSR, RHS, Solution, iRHS, iSolution, NULL, NULL,
(int)Matrix->SMPkluMatrix->KLUmatrixN, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
} else {
Ax_CSR = (double *) malloc ((size_t)Matrix->SMPkluMatrix->KLUmatrixNZ * sizeof (double)) ;
klu_convert_matrix_in_CSR (Matrix->SMPkluMatrix->KLUmatrixAp, Matrix->SMPkluMatrix->KLUmatrixAi, Matrix->SMPkluMatrix->KLUmatrixAx, Ap_CSR, Ai_CSR,
Ax_CSR, (int)Matrix->SMPkluMatrix->KLUmatrixN, (int)Matrix->SMPkluMatrix->KLUmatrixNZ, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
klu_matrix_vector_multiply (Ap_CSR, Ai_CSR, Ax_CSR, RHS, Solution, NULL, NULL,
(int)Matrix->SMPkluMatrix->KLUmatrixN, Matrix->SMPkluMatrix->KLUmatrixCommon) ;
iSolution = iRHS ;
}
free (Ap_CSR) ;
free (Ai_CSR) ;
free (Ax_CSR) ;
} else {
spMultiply (Matrix->SPmatrix, RHS, Solution, iRHS, iSolution) ;
}
}
Reference in New Issue View Git Blame Copy Permalink