/* README :
* The routines in this file use a pointer to a Matrix
* at the Fortran Level this pointer is stored as an integer
* so we cast pointer to long int
* f_(fmat)
* long *fmat
* *fmat = (long)spCreate(*n,0,&error);
*
* At fortran level the integer must be an integer *4
* in order to store a C long
* since we are in scilab
* the pointer transmitted to f_ is an istk(il1) it can in fact contain
* something as long as a double
* Copyright ENPC (Chancelier)
*/
/*
* IMPORTS
*
* >>> Import descriptions:
* spConfig.h
* Macros that customize the sparse matrix routines.
* spmatrix.h
* Macros and declarations to be imported by the user.
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#define spINSIDE_SPARSE
#include "spConfig.h"
#include "spmatrix.h"
#include "spDefs.h"
#include "../machine.h"
/*
*
* lufact1 >>> Creation and LU factorisation of a sparse matrix
* Entry <<
* val,rc
* arrays of size k and kx2 a(rc(i),rc(*k+i))=val(i)
* for i=0,(*k-1)
* *n : size of the created square matrix
* *k : number of given values
* *eps : The machine precision number
* Return >>
* fmat : pointer to a long int which is a cast of an adress
* *nrank : The numerical rank
*/
static void
spFixThresold(eMatrix,eps,releps)
char *eMatrix;
double eps,releps;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
Matrix->AbsThreshold = eps;
Matrix->RelThreshold = releps;
}
static void
spGetNumRank(eMatrix,n)
char *eMatrix;
int *n;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
*n = Matrix->NumRank;
}
extern void cerro();
void
C2F(lufact1)(val,lln,col,n,nel,fmat,eps,releps,nrank,ierr)
double *val,*eps,*releps;
long *fmat;
int *n,*nel,*nrank,*lln,*col,*ierr;
{
int error,i,i0,i1,k,j;
spREAL *pelement;
*ierr = 0;
*fmat = (long)spCreate(*n,0,&error);
if (error != spOKAY) {
*ierr = 1;
return;
}
i0=0;
i1=i0;
i=1;
for (k = 0 ;k < *nel; k++) {
i0=i0+1;
while (i0-i1 > lln[i-1]) {
i1=i0;
i=i+1;
i0=i0+1;
}
j=col[k];
pelement = spGetElement((char*) *fmat,i,j);
if (pelement == 0) {
*ierr=2;
return;
}
spADD_REAL_ELEMENT(pelement,(spREAL)(val[k]));
}
/* Fix the AbsThresold with scilex %eps */
spFixThresold((char*) *fmat,*eps,*releps);
/* spPrint((char *) *fmat,1,1,1); */
error = spFactor((char*) *fmat);
spGetNumRank((char *) *fmat,nrank);
switch (error) {
case spZERO_DIAG:
cerro("zero_diag: A zero was encountered on the diagonal the matrix ");
break;
case spNO_MEMORY:
*ierr=3;
break;
case spSINGULAR:
*ierr=-1; /*Singular matrix" */
break;
case spSMALL_PIVOT:
*ierr=-2; /* matrix is singular at precision level */
break;
}
}
/*
* lusolve1 >>> Solves fmat*x=b
* *fmat : a pointer to the sparse matrix factored by lufact
* b,v
* two arrays of size n the matrix size
*/
extern void Cout(char *str);
static int
spSolveCheck(eMatrix)
char *eMatrix;
{
MatrixPtr Matrix = (MatrixPtr) eMatrix;
if (Matrix->Error == spSINGULAR OR Matrix->Error == spSMALL_PIVOT )
{
Cout("A matrix is singular or ill conditioned ");
Cout("result will be fine only if b is in Im(A)");
}
return(1);
}
void C2F(lusolve1)(fmat,b,x)
double *b, *x;
long *fmat;
{
/* if (spSolveCheck((char *) *fmat) != 0)*/
spSolve((char*) *fmat,(spREAL*)b,(spREAL*)x);
}
/*
* ludel1 >>> delete sparse matrix
* *fmat : a pointer to the sparse matrix factored by lufact
*/
void C2F(ludel1)(fmat)
long *fmat;
{
spDestroy((char*) *fmat);
}
/*
* lusize >>> returns in n the size of the sparse matrix
* *fmat : a pointer to the sparse matrix factored by lufact
*/
static void
spSize(eMatrix,n)
char *eMatrix;
int *n;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
*n=Matrix->Size;
}
void C2F(lusize)(fmat,n)
long *fmat;
int *n;
{
spSize((char *) *fmat,n);
}
/*
* luget1 >>> extract the LU coded matrix into a full array
* sigg,sigd :
* two arrays of size n which code permutations
* lu :
* an array coded matrix of size nxn where lu will be stored
*/
/* filling right permutation */
GetSigD(Matrix,indsigd,sigd)
MatrixPtr Matrix;
double sigd[];
int indsigd[];
{
int I,J,mc=0,last=0;
int Size=Matrix->Size;
for (I = 1; I <= Size; I++)
{
indsigd[I-1]=1;
indsigd[Size+I-1]= Matrix->IntToExtColMap[I];
sigd[I-1]=1.0;
}
/* counting missing colums */
for (I = 1; I <= Size; I++)
if (Matrix->ExtToIntColMap[I]== -1) mc++;
/* filling missing colums */
if (mc != 0)
{
for (I = Size -(mc)+1 ; I <= Size; I++)
{
for ( J=last+1; J <=Size; J++)
{
if (Matrix->ExtToIntColMap[J]==-1)
{
last=J;break;
}
}
indsigd[I-1]=1;
indsigd[Size+I-1]= last;
}
}
}
/* filling left permutation */
GetSigG(Matrix,indsigg,sigg)
MatrixPtr Matrix;
double sigg[];
int indsigg[];
{
int Size=Matrix->Size;
int I,J,mc=0,last=0;
/* counting missing Rows*/
for (I = 1; I <= Size; I++)
if (Matrix->ExtToIntRowMap[I]== -1) mc++;
for (I = 1; I <= Size-mc ; I++)
{
indsigg[I-1]= 1;
indsigg[Size+Matrix->IntToExtRowMap[I]-1]=I;
sigg[I-1]=1;
}
/* filling missing Rows */
if (mc != 0)
{
for (I = Size -(mc)+1 ; I <= Size; I++)
{
for ( J=last+1; J <=Size; J++)
{
if (Matrix->ExtToIntRowMap[J]==-1)
{
last=J;break;
}
}
indsigg[I-1]= 1;
indsigg[Size+last-1]=I;
sigg[I-1]=1;
}
}
}
static void
spLuget(eMatrix,indP,P,indl,l,indu,u,indQ,Q)
char *eMatrix;
int *indP,*indl,*indu,*indQ;
double *P,*Q,*l,*u;
{
int I,J;
int lsize,usize;
MatrixPtr Matrix = (MatrixPtr) eMatrix;
ElementPtr pElement;
int Size;
Size = Matrix->Size;
GetSigD(Matrix,indQ,Q);
GetSigG(Matrix,indP,P);
for (J = 1; J <= Size ; J++)
{
indl[J-1] = 0;
indu[J-1] = 0;
}
lsize=0;
usize=0;
for (I = 1; I <= Size ; I++)
{
indu[I-1]=indu[I-1]+1;
indu[Size+usize]=I;
u[usize]=1.0;
usize=usize+1;
pElement = Matrix->FirstInRow[I];
while ( pElement != NULL )
{
J = pElement->Col;
if (I >= J) {
indl[I-1] = indl[I-1]+1;
indl[Size+lsize]=J;
l[lsize]=(double) pElement->Real ;
lsize=lsize+1;
}
else {
indu[I-1] = indu[I-1]+1;
indu[Size+usize]=J;
u[usize]=(double) pElement->Real ;
usize=usize+1;
}
pElement = pElement->NextInRow;
};
};
}
void C2F(luget1)(fmat,indP,P,indl,l,indu,u,indQ,Q)
long *fmat;
double *P,*Q,*l,*u;
int *indP,*indl,*indu,*indQ;
{
spLuget((char *) *fmat,indP,P,indl,l,indu,u,indQ,Q);
}
/*
* lusiz1 >>> extract the L and U number of non zero elements
* lsize and usize
*/
static void
spLusiz(eMatrix,lsize,usize)
char *eMatrix;
int *lsize,*usize;
{
int J;
MatrixPtr Matrix = (MatrixPtr) eMatrix;
ElementPtr pElement;
int Size;
Size = Matrix->Size;
*lsize=0;
*usize=Size;
for (J = 1; J <= Size ; J++)
{
pElement = Matrix->FirstInCol[J];
while ( pElement != NULL )
{
if (pElement->Row >= J)
*lsize=*lsize+1;
else
*usize=*usize+1;
pElement = pElement->NextInCol;
};
};
}
void C2F(lusiz1)(fmat,lsize,usize)
long *fmat;
int *lsize,*usize;
{
spLusiz((char *) *fmat,lsize,usize);
}
syntax highlighted by Code2HTML, v. 0.9.1