# 1  "spfactor.c"
static char copyright[] = "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert";
static char RCSid[] = "@(#)$Header: spFactor.c,v 1.2 88/06/18 11:14:10 kundert Exp $";
struct spTemplate {
    double  *Element1;
    double  *Element2;
    double  *Element3Negated;
    double  *Element4Negated;
    } ;
extern void spClear();
extern double spCondition();
extern char  *spCreate();
extern void spDeleteRowAndCol();
extern void spDestroy();
extern void spDeterminant();
extern int spElementCount();
extern int spError();
extern int spFactor();
extern int spFileMatrix();
extern int spFileStats();
extern int spFileVector();
extern int spFillinCount();
extern int spGetAdmittance();
extern double  *spGetElement();
extern char  *spGetInitInfo();
extern int spGetOnes();
extern int spGetQuad();
extern int spGetSize();
extern int spInitialize();
extern void spInstallInitInfo();
extern double spLargestElement();
extern void spMNA_Preorder();
extern void spMultiply();
extern void spMultTransposed();
extern double spNorm();
extern int spOrderAndFactor();
extern void spPartition();
extern void spPrint();
extern double spPseudoCondition();
extern double spRoundoff();
extern void spScale();
extern void spSetComplex();
extern void spSetReal();
extern void spSolve();
extern void spSolveTransposed();
extern void spStripFills();
extern void spWhereSingular();
typedef struct {
    int _cnt;
    unsigned char  *_ptr;
    unsigned char  *_base;
    char _flag;
    char _file;
    }  FILE;
extern struct usptr_s  *_sproced;
extern FILE _iob[100];
extern FILE  *fopen();
extern FILE  *fdopen();
extern FILE  *freopen();
extern FILE  *popen();
extern FILE  *tmpfile();
extern long ftell();
extern void rewind();
extern void setbuf();
extern char  *ctermid();
extern char  *cuserid();
extern char  *fgets();
extern char  *gets();
extern char  *tempnam();
extern char  *tmpnam();
extern int fclose();
extern int fflush();
extern int fread();
extern int fwrite();
extern int fseek();
extern int fgetc();
extern int getw();
extern int pclose();
extern int printf();
extern int fprintf();
extern int sprintf();
extern int vprintf();
extern int vfprintf();
extern int vsprintf();
extern int fputc();
extern int putw();
extern int puts();
extern int fputs();
extern int scanf();
extern int fscanf();
extern int sscanf();
extern int setvbuf();
extern int system();
extern int ungetc();
extern unsigned char  *_bufendtab[];
extern char  *malloc();
extern char  *calloc();
extern char  *realloc();
extern int free();
extern int abort();
typedef double RealNumber;
typedef double  *RealVector;
typedef struct {
    double Real;
    double Imag;
    }  ComplexNumber,  *ComplexVector;
struct MatrixElement {
    double Real;
    double Imag;
    int Row;
    int Col;
    struct MatrixElement  *NextInRow;
    struct MatrixElement  *NextInCol;
    char  *pInitInfo;
    } ;
typedef struct MatrixElement  *ElementPtr;
typedef ElementPtr  *ArrayOfElementPtrs;
struct AllocationRecord {
    char  *AllocatedPtr;
    struct AllocationRecord  *NextRecord;
    } ;
typedef struct AllocationRecord  *AllocationListPtr;
struct FillinListNodeStruct {
    ElementPtr pFillinList;
    int NumberOfFillinsInList;
    struct FillinListNodeStruct  *Next;
    } ;
struct MatrixFrame {
    double AbsThreshold;
    int AllocatedSize;
    int AllocatedExtSize;
    int Complex;
    int CurrentSize;
    ArrayOfElementPtrs Diag;
    int  *DoCmplxDirect;
    int  *DoRealDirect;
    int Elements;
    int Error;
    int ExtSize;
    int  *ExtToIntColMap;
    int  *ExtToIntRowMap;
    int Factored;
    int Fillins;
    ArrayOfElementPtrs FirstInCol;
    ArrayOfElementPtrs FirstInRow;
    unsigned long ID;
    RealVector Intermediate;
    int InternalVectorsAllocated;
    int  *IntToExtColMap;
    int  *IntToExtRowMap;
    int  *MarkowitzRow;
    int  *MarkowitzCol;
    long  *MarkowitzProd;
    int MaxRowCountInLowerTri;
    int NeedsOrdering;
    int NumberOfInterchangesIsOdd;
    int Partitioned;
    int PivotsOriginalCol;
    int PivotsOriginalRow;
    char PivotSelectionMethod;
    int PreviousMatrixWasComplex;
    double RelThreshold;
    int Reordered;
    int RowsLinked;
    int SingularCol;
    int SingularRow;
    int Singletons;
    int Size;
    struct MatrixElement TrashCan;
    AllocationListPtr TopOfAllocationList;
    int RecordsRemaining;
    ElementPtr NextAvailElement;
    int ElementsRemaining;
    ElementPtr NextAvailFillin;
    int FillinsRemaining;
    struct FillinListNodeStruct  *FirstFillinListNode;
    struct FillinListNodeStruct  *LastFillinListNode;
    } ;
typedef struct MatrixFrame  *MatrixPtr;
# 173  "spfactor.c"
int spOrderAndFactor( eMatrix, RHS, RelThreshold, AbsThreshold, DiagPivoting )
char  *eMatrix;
    double  *RHS;
    double RelThreshold;
    double AbsThreshold;
    int DiagPivoting;
 
{
    MatrixPtr Matrix = (MatrixPtr )(eMatrix);
    ElementPtr pPivot;
    int Step;
    int Size;
    int ReorderingRequired;
    ElementPtr SearchForPivot();
    double LargestInCol;
    double FindLargestInCol();
    register double DD1;
    register double DD2;
 
# 186  "spfactor.c"
    if (!((Matrix != 0 && (*Matrix).ID == 0x772773 && (*Matrix).Error >= 0 && (*Matrix).Error < 2) && !(*Matrix).Factored)) {
# 186  "spfactor.c"
        fflush( (&_iob[1]) );
# 186  "spfactor.c"
        fprintf( (&_iob[2]), "sparse: panic in file `%s' at line %d.\n", "spfactor.c", 186 );
# 186  "spfactor.c"
        fflush( (&_iob[2]) );
# 186  "spfactor.c"
        abort(  );
# 186  "spfactor.c"
        }
# 186  "spfactor.c"
    ;
 
    (*Matrix).Error = 0;
    Size = (*Matrix).Size;
    if (RelThreshold <= 0.0) {
# 190  "spfactor.c"
        RelThreshold = (*Matrix).RelThreshold;
# 190  "spfactor.c"
        }
    if (RelThreshold > 1.0) {
# 191  "spfactor.c"
        RelThreshold = (*Matrix).RelThreshold;
# 191  "spfactor.c"
        }
    (*Matrix).RelThreshold = RelThreshold;
    if (AbsThreshold < 0.0) {
# 193  "spfactor.c"
        AbsThreshold = (*Matrix).AbsThreshold;
# 193  "spfactor.c"
        }
    (*Matrix).AbsThreshold = AbsThreshold;
    ReorderingRequired = 0;
 
    if (!(*Matrix).NeedsOrdering) {
# 200  "spfactor.c"
        Step = 1;
# 200  "spfactor.c"
        while ( Step <= Size ) {
            pPivot = (*Matrix).Diag[Step];
            LargestInCol = FindLargestInCol( (*pPivot).NextInCol );
            if (((*pPivot).Real) < 0.0) {
# 203  "spfactor.c"
                DD1 = -((*pPivot).Real);
# 203  "spfactor.c"
            } else {
# 203  "spfactor.c"
                DD1 = ((*pPivot).Real);
# 203  "spfactor.c"
                }
# 203  "spfactor.c"
            if (((*pPivot).Imag) < 0.0) {
# 203  "spfactor.c"
                DD2 = -((*pPivot).Imag);
# 203  "spfactor.c"
            } else {
# 203  "spfactor.c"
                DD2 = ((*pPivot).Imag);
# 203  "spfactor.c"
                }
# 203  "spfactor.c"
            if ((LargestInCol * RelThreshold < (DD1 + DD2))) {
                if ((*Matrix).Complex) {
                    ComplexRowColElimination( Matrix, pPivot );
# 205  "spfactor.c"
                } else {
 
                    RealRowColElimination( Matrix, pPivot );
# 207  "spfactor.c"
                    }
# 207  "spfactor.c"
            } else {
# 210  "spfactor.c"
                ReorderingRequired = 1;
                goto l2;
# 211  "spfactor.c"
                }
# 200  "spfactor.c"
            Step++;
# 200  "spfactor.c"
            }
# 211  "spfactor.c"
l2:
# 214  "spfactor.c"
        if (!ReorderingRequired) {
            goto Done;
# 215  "spfactor.c"
            }
# 215  "spfactor.c"
        ;
# 215  "spfactor.c"
    } else {
# 236  "spfactor.c"
        Step = 1;
        if (!(*Matrix).RowsLinked) {
            spcLinkRows( Matrix );
# 238  "spfactor.c"
            }
        if (!(*Matrix).InternalVectorsAllocated) {
            CreateInternalVectors( Matrix );
# 240  "spfactor.c"
            }
        if ((*Matrix).Error >= 2) {
            return (*Matrix).Error;
# 242  "spfactor.c"
            }
# 242  "spfactor.c"
        }
# 246  "spfactor.c"
    CountMarkowitz( Matrix, RHS, Step );
    MarkowitzProducts( Matrix, Step );
    (*Matrix).MaxRowCountInLowerTri = -1;
# 251  "spfactor.c"
    for ( ; Step<=Size; Step++ ) {
        pPivot = SearchForPivot( Matrix, Step, DiagPivoting );
        if (pPivot == 0) {
# 253  "spfactor.c"
            return MatrixIsSingular( Matrix, Step );
# 253  "spfactor.c"
            }
        ExchangeRowsAndCols( Matrix, pPivot, Step );
 
        if ((*Matrix).Complex) {
            ComplexRowColElimination( Matrix, pPivot );
# 257  "spfactor.c"
        } else {
 
            RealRowColElimination( Matrix, pPivot );
# 259  "spfactor.c"
            }
 
        if ((*Matrix).Error >= 2) {
# 261  "spfactor.c"
            return (*Matrix).Error;
# 261  "spfactor.c"
            }
        UpdateMarkowitzNumbers( Matrix, pPivot );
    }
# 269  "spfactor.c"
Done:
    (*Matrix).NeedsOrdering = 0;
    (*Matrix).Reordered = 1;
    (*Matrix).Factored = 1;
 
    return (*Matrix).Error;
}
# 1  "spfactor.c"
# 315  "spfactor.c"
int spFactor( eMatrix )
char  *eMatrix;
 
{
    MatrixPtr Matrix = (MatrixPtr )(eMatrix);
    register ElementPtr pElement;
    register ElementPtr pColumn;
    register int Step;
    register int Size;
    double Mult;
    register int II1;
    register double  *Dest;
    register double  * *pDest;
 
# 326  "spfactor.c"
    if (!((Matrix != 0 && (*Matrix).ID == 0x772773 && (*Matrix).Error >= 0 && (*Matrix).Error < 2) && !(*Matrix).Factored)) {
# 326  "spfactor.c"
        fflush( (&_iob[1]) );
# 326  "spfactor.c"
        fprintf( (&_iob[2]), "sparse: panic in file `%s' at line %d.\n", "spfactor.c", 326 );
# 326  "spfactor.c"
        fflush( (&_iob[2]) );
# 326  "spfactor.c"
        abort(  );
# 326  "spfactor.c"
        }
# 326  "spfactor.c"
    ;
 
    if ((*Matrix).NeedsOrdering) {
        return spOrderAndFactor( eMatrix, (RealVector )(0), 0.0, 0.0, 1 );
# 329  "spfactor.c"
        }
# 332  "spfactor.c"
    if (!(*Matrix).Partitioned) {
# 332  "spfactor.c"
        spPartition( eMatrix, 0 );
# 332  "spfactor.c"
        }
 
    if ((*Matrix).Complex) {
# 334  "spfactor.c"
        return FactorComplexMatrix( Matrix );
# 334  "spfactor.c"
        }
# 338  "spfactor.c"
    Size = (*Matrix).Size;
 
    if ((*(*Matrix).Diag[1]).Real == 0.0) {
# 340  "spfactor.c"
        return ZeroPivot( Matrix, 1 );
# 340  "spfactor.c"
        }
    (*(*Matrix).Diag[1]).Real = 1.0 / (*(*Matrix).Diag[1]).Real;
# 344  "spfactor.c"
    for ( Step = 2; Step<=Size; Step++ ) {
        if ((*Matrix).DoRealDirect[Step]) {
# 375  "spfactor.c"
            Dest = (double  *)((*Matrix).Intermediate);
# 350  "spfactor.c"
            pElement = (*Matrix).FirstInCol[Step];
            while ( pElement != 0 ) {
                Dest[(*pElement).Row] = (*pElement).Real;
                pElement = (*pElement).NextInCol;
# 353  "spfactor.c"
                }
# 357  "spfactor.c"
            pColumn = (*Matrix).FirstInCol[Step];
            while ( (*pColumn).Row < Step ) {
                pElement = (*Matrix).Diag[(*pColumn).Row];
                (*pColumn).Real = Dest[(*pColumn).Row] * (*pElement).Real;
                pElement = (*pElement).NextInCol;
# 361  "spfactor.c"
                while ( pElement != 0 ) {
                    Dest[(*pElement).Row] -=  (*pColumn).Real * (*pElement).Real;
# 361  "spfactor.c"
                    pElement = (*pElement).NextInCol;
# 361  "spfactor.c"
                    }
 
                pColumn = (*pColumn).NextInCol;
# 363  "spfactor.c"
                }
# 367  "spfactor.c"
            pElement = (*(*Matrix).Diag[Step]).NextInCol;
            while ( pElement != 0 ) {
                (*pElement).Real = Dest[(*pElement).Row];
                pElement = (*pElement).NextInCol;
# 370  "spfactor.c"
                }
# 374  "spfactor.c"
            if (Dest[Step] == 0.0) {
# 374  "spfactor.c"
                return ZeroPivot( Matrix, Step );
# 374  "spfactor.c"
                }
            (*(*Matrix).Diag[Step]).Real = 1.0 / Dest[Step];
# 375  "spfactor.c"
        } else {
# 401  "spfactor.c"
            pDest = (double  * *)((*Matrix).Intermediate);
# 382  "spfactor.c"
            pElement = (*Matrix).FirstInCol[Step];
            while ( pElement != 0 ) {
                pDest[(*pElement).Row] = &(*pElement).Real;
                pElement = (*pElement).NextInCol;
# 385  "spfactor.c"
                }
# 389  "spfactor.c"
            pColumn = (*Matrix).FirstInCol[Step];
            while ( (*pColumn).Row < Step ) {
                pElement = (*Matrix).Diag[(*pColumn).Row];
                *pDest[(*pColumn).Row] *=  (*pElement).Real;
# 392  "spfactor.c"
                Mult = (*pDest[(*pColumn).Row]);
                pElement = (*pElement).NextInCol;
# 393  "spfactor.c"
                while ( pElement != 0 ) {
                    *pDest[(*pElement).Row] -=  Mult * (*pElement).Real;
# 393  "spfactor.c"
                    pElement = (*pElement).NextInCol;
# 393  "spfactor.c"
                    }
 
                pColumn = (*pColumn).NextInCol;
# 395  "spfactor.c"
                }
# 399  "spfactor.c"
            if ((*(*Matrix).Diag[Step]).Real == 0.0) {
                return ZeroPivot( Matrix, Step );
# 400  "spfactor.c"
                }
            (*(*Matrix).Diag[Step]).Real = 1.0 / (*(*Matrix).Diag[Step]).Real;
# 401  "spfactor.c"
            }
    }
# 405  "spfactor.c"
    (*Matrix).Factored = 1;
    II1 = 0;
# 406  "spfactor.c"
    (*Matrix).Error = II1;
# 406  "spfactor.c"
    return II1;
}
# 1  "spfactor.c"
# 435  "spfactor.c"
static int FactorComplexMatrix( Matrix )
MatrixPtr Matrix;
 
{
    register ElementPtr pElement;
    register ElementPtr pColumn;
    register int Step;
    register int Size;
    ComplexNumber Mult;
    ComplexNumber Pivot;
    register double DD1;
    register double DD2;
    register int II1;
    double DD3;
    double DD6;
    register ComplexNumber  *Dest;
    register double DD4;
    register double DD5;
    double r_;
    register ComplexNumber  * *pDest;
    register double DD7;
    register double DD8;
 
# 445  "spfactor.c"
    if (!((*Matrix).Complex)) {
# 445  "spfactor.c"
        fflush( (&_iob[1]) );
# 445  "spfactor.c"
        fprintf( (&_iob[2]), "sparse: panic in file `%s' at line %d.\n", "spfactor.c", 445 );
# 445  "spfactor.c"
        fflush( (&_iob[2]) );
# 445  "spfactor.c"
        abort(  );
# 445  "spfactor.c"
        }
# 445  "spfactor.c"
    ;
 
    Size = (*Matrix).Size;
    pElement = (*Matrix).Diag[1];
    if (((*pElement).Real) < 0.0) {
# 449  "spfactor.c"
        DD1 = -((*pElement).Real);
# 449  "spfactor.c"
    } else {
# 449  "spfactor.c"
        DD1 = ((*pElement).Real);
# 449  "spfactor.c"
        }
# 449  "spfactor.c"
    if (((*pElement).Imag) < 0.0) {
# 449  "spfactor.c"
        DD2 = -((*pElement).Imag);
# 449  "spfactor.c"
    } else {
# 449  "spfactor.c"
        DD2 = ((*pElement).Imag);
# 449  "spfactor.c"
        }
# 449  "spfactor.c"
    if ((DD1 + DD2) == 0.0) {
# 449  "spfactor.c"
        return ZeroPivot( Matrix, 1 );
# 449  "spfactor.c"
        }
 
    if (((*pElement).Real >= (*pElement).Imag && (*pElement).Real > -(*pElement).Imag) || ((*pElement).Real < (*pElement).Imag &&
        (*pElement).Real <= -(*pElement).Imag)) {
# 451  "spfactor.c"
        DD3 = (*pElement).Imag / (*pElement).Real;
# 451  "spfactor.c"
        (*pElement).Real = 1.0 / ((*pElement).Real + DD3 * (*pElement).Imag);
# 451  "spfactor.c"
        (*pElement).Imag = -DD3 * ((*pElement).Real);
# 451  "spfactor.c"
    } else {
# 451  "spfactor.c"
        DD3 = (*pElement).Real / (*pElement).Imag;
# 451  "spfactor.c"
        (*pElement).Imag = -1.0 / ((*pElement).Imag + DD3 * (*pElement).Real);
# 451  "spfactor.c"
        (*pElement).Real = -DD3 * ((*pElement).Imag);
# 451  "spfactor.c"
        }
# 451  "spfactor.c"
    ;
# 454  "spfactor.c"
    for ( Step = 2; Step<=Size; Step++ ) {
        if ((*Matrix).DoCmplxDirect[Step]) {
# 458  "spfactor.c"
            Dest = (ComplexNumber  *)((*Matrix).Intermediate);
# 461  "spfactor.c"
            pElement = (*Matrix).FirstInCol[Step];
            while ( pElement != 0 ) {
                Dest[(*pElement).Row] = *(ComplexNumber  *)(pElement);
                pElement = (*pElement).NextInCol;
# 464  "spfactor.c"
                }
# 468  "spfactor.c"
            pColumn = (*Matrix).FirstInCol[Step];
            while ( (*pColumn).Row < Step ) {
                pElement = (*Matrix).Diag[(*pColumn).Row];
 
                Mult.Real = Dest[(*pColumn).Row].Real * (*pElement).Real - Dest[(*pColumn).Row].Imag * (*pElement).Imag;
# 472  "spfactor.c"
                Mult.Imag = Dest[(*pColumn).Row].Real * (*pElement).Imag + Dest[(*pColumn).Row].Imag * (*pElement).Real;
# 472  "spfactor.c"
                ;
                (*pColumn).Real = Mult.Real;
# 473  "spfactor.c"
                (*pColumn).Imag = Mult.Imag;
# 473  "spfactor.c"
                ;
                pElement = (*pElement).NextInCol;
# 474  "spfactor.c"
                while ( pElement != 0 ) {
 
                    Dest[(*pElement).Row].Real -=  Mult.Real * (*pElement).Real - Mult.Imag * (*pElement).Imag;
# 476  "spfactor.c"
                    Dest[(*pElement).Row].Imag -=  Mult.Real * (*pElement).Imag + Mult.Imag * (*pElement).Real;
# 476  "spfactor.c"
                    ;
# 474  "spfactor.c"
                    pElement = (*pElement).NextInCol;
# 474  "spfactor.c"
                    }
# 478  "spfactor.c"
                pColumn = (*pColumn).NextInCol;
# 478  "spfactor.c"
                }
# 482  "spfactor.c"
            pElement = (*(*Matrix).Diag[Step]).NextInCol;
            while ( pElement != 0 ) {
                *(ComplexNumber  *)(pElement) = Dest[(*pElement).Row];
                pElement = (*pElement).NextInCol;
# 485  "spfactor.c"
                }
# 489  "spfactor.c"
            Pivot = Dest[Step];
            if ((Pivot.Real) < 0.0) {
# 490  "spfactor.c"
                DD4 = -(Pivot.Real);
# 490  "spfactor.c"
            } else {
# 490  "spfactor.c"
                DD4 = (Pivot.Real);
# 490  "spfactor.c"
                }
# 490  "spfactor.c"
            if ((Pivot.Imag) < 0.0) {
# 490  "spfactor.c"
                DD5 = -(Pivot.Imag);
# 490  "spfactor.c"
            } else {
# 490  "spfactor.c"
                DD5 = (Pivot.Imag);
# 490  "spfactor.c"
                }
# 490  "spfactor.c"
            if ((DD4 + DD5) == 0.0) {
# 490  "spfactor.c"
                return ZeroPivot( Matrix, Step );
# 490  "spfactor.c"
                }
            if ((Pivot.Real >= Pivot.Imag && Pivot.Real > -Pivot.Imag) || (Pivot.Real < Pivot.Imag && Pivot.Real <= -Pivot.Imag)
                ) {
# 491  "spfactor.c"
                DD6 = Pivot.Imag / Pivot.Real;
# 491  "spfactor.c"
                (*(*Matrix).Diag[Step]).Real = 1.0 / (Pivot.Real + DD6 * Pivot.Imag);
# 491  "spfactor.c"
                (*(*Matrix).Diag[Step]).Imag = -DD6 * ((*(*Matrix).Diag[Step]).Real);
# 491  "spfactor.c"
            } else {
# 491  "spfactor.c"
                DD6 = Pivot.Real / Pivot.Imag;
# 491  "spfactor.c"
                (*(*Matrix).Diag[Step]).Imag = -1.0 / (Pivot.Imag + DD6 * Pivot.Real);
# 491  "spfactor.c"
                (*(*Matrix).Diag[Step]).Real = -DD6 * ((*(*Matrix).Diag[Step]).Imag);
# 491  "spfactor.c"
                }
# 491  "spfactor.c"
            ;
# 491  "spfactor.c"
        } else {
# 496  "spfactor.c"
            pDest = (ComplexNumber  * *)((*Matrix).Intermediate);
# 499  "spfactor.c"
            pElement = (*Matrix).FirstInCol[Step];
            while ( pElement != 0 ) {
                pDest[(*pElement).Row] = (ComplexNumber  *)(pElement);
                pElement = (*pElement).NextInCol;
# 502  "spfactor.c"
                }
# 506  "spfactor.c"
            pColumn = (*Matrix).FirstInCol[Step];
            while ( (*pColumn).Row < Step ) {
                pElement = (*Matrix).Diag[(*pColumn).Row];
 
                Mult.Real = (*pDest[(*pColumn).Row]).Real * (*pElement).Real - (*pDest[(*pColumn).Row]).Imag * (*pElement).Imag;
# 510  "spfactor.c"
                Mult.Imag = (*pDest[(*pColumn).Row]).Real * (*pElement).Imag + (*pDest[(*pColumn).Row]).Imag * (*pElement).Real;
# 510  "spfactor.c"
                ;
                (*pDest[(*pColumn).Row]).Real = Mult.Real;
# 511  "spfactor.c"
                (*pDest[(*pColumn).Row]).Imag = Mult.Imag;
# 511  "spfactor.c"
                ;
                pElement = (*pElement).NextInCol;
# 512  "spfactor.c"
                while ( pElement != 0 ) {
 
                    (*pDest[(*pElement).Row]).Real -=  Mult.Real * (*pElement).Real - Mult.Imag * (*pElement).Imag;
# 514  "spfactor.c"
                    (*pDest[(*pElement).Row]).Imag -=  Mult.Real * (*pElement).Imag + Mult.Imag * (*pElement).Real;
# 514  "spfactor.c"
                    ;
# 512  "spfactor.c"
                    pElement = (*pElement).NextInCol;
# 512  "spfactor.c"
                    }
# 516  "spfactor.c"
                pColumn = (*pColumn).NextInCol;
# 516  "spfactor.c"
                }
# 520  "spfactor.c"
            pElement = (*Matrix).Diag[Step];
            if (((*pElement).Real) < 0.0) {
# 521  "spfactor.c"
                DD7 = -((*pElement).Real);
# 521  "spfactor.c"
            } else {
# 521  "spfactor.c"
                DD7 = ((*pElement).Real);
# 521  "spfactor.c"
                }
# 521  "spfactor.c"
            if (((*pElement).Imag) < 0.0) {
# 521  "spfactor.c"
                DD8 = -((*pElement).Imag);
# 521  "spfactor.c"
            } else {
# 521  "spfactor.c"
                DD8 = ((*pElement).Imag);
# 521  "spfactor.c"
                }
# 521  "spfactor.c"
            if ((DD7 + DD8) == 0.0) {
# 521  "spfactor.c"
                return ZeroPivot( Matrix, Step );
# 521  "spfactor.c"
                }
            if (((*pElement).Real >= (*pElement).Imag && (*pElement).Real > -(*pElement).Imag) || ((*pElement).Real < (*pElement).
                Imag && (*pElement).Real <= -(*pElement).Imag)) {
# 522  "spfactor.c"
                r_ = (*pElement).Imag / (*pElement).Real;
# 522  "spfactor.c"
                (*pElement).Real = 1.0 / ((*pElement).Real + r_ * (*pElement).Imag);
# 522  "spfactor.c"
                (*pElement).Imag = -r_ * ((*pElement).Real);
# 522  "spfactor.c"
            } else {
# 522  "spfactor.c"
                r_ = (*pElement).Real / (*pElement).Imag;
# 522  "spfactor.c"
                (*pElement).Imag = -1.0 / ((*pElement).Imag + r_ * (*pElement).Real);
# 522  "spfactor.c"
                (*pElement).Real = -r_ * ((*pElement).Imag);
# 522  "spfactor.c"
                }
# 522  "spfactor.c"
            ;
# 522  "spfactor.c"
            }
    }
# 526  "spfactor.c"
    (*Matrix).Factored = 1;
    II1 = 0;
# 527  "spfactor.c"
    (*Matrix).Error = II1;
# 527  "spfactor.c"
    return II1;
}
# 1  "spfactor.c"
# 574  "spfactor.c"
void spPartition( eMatrix, Mode )
char  *eMatrix;
    int Mode;
 
{
    MatrixPtr Matrix = (MatrixPtr )(eMatrix);
    register ElementPtr pElement;
    register ElementPtr pColumn;
    register int Step;
    register int Size;
    register int  *Nc;
    register int  *No;
    register int  *Nm;
    int  *DoRealDirect;
    int  *DoCmplxDirect;
    int Step1;
    int Step2;
 
# 586  "spfactor.c"
    if (!(Matrix != 0 && (*Matrix).ID == 0x772773)) {
# 586  "spfactor.c"
        fflush( (&_iob[1]) );
# 586  "spfactor.c"
        fprintf( (&_iob[2]), "sparse: panic in file `%s' at line %d.\n", "spfactor.c", 586 );
# 586  "spfactor.c"
        fflush( (&_iob[2]) );
# 586  "spfactor.c"
        abort(  );
# 586  "spfactor.c"
        }
# 586  "spfactor.c"
    ;
    if ((*Matrix).Partitioned) {
# 587  "spfactor.c"
        return ;
# 587  "spfactor.c"
        }
    Size = (*Matrix).Size;
    DoRealDirect = (*Matrix).DoRealDirect;
    DoCmplxDirect = (*Matrix).DoCmplxDirect;
    (*Matrix).Partitioned = 1;
# 594  "spfactor.c"
    if (Mode == 0) {
# 594  "spfactor.c"
        Mode = 3;
# 594  "spfactor.c"
        }
    if (Mode == 1) {
#pragma parallel if(Size > 333) byvalue(Size) shared(Step, DoRealDirect) local(Step1)
        {
# 596  "spfactor.c"
#pragma pfor iterate(Step1=1;Size;1)
# 596  "spfactor.c"
            for ( Step1 = 1; Step1<=Size; Step1++ ) {
 
                DoRealDirect[Step1] = 1;
            }
# 598  "spfactor.c"
        }
# 596  "spfactor.c"
        Step = 1 + ((Size)>(0) ? (Size) : (0));
# 601  "spfactor.c"
        DoCmplxDirect[Step] = 1;
 
        return ;
# 603  "spfactor.c"
        }
 
    if (Mode == 2) {
#pragma parallel if(Size > 333) byvalue(Size) shared(Step, DoRealDirect) local(Step2)
        {
# 606  "spfactor.c"
#pragma pfor iterate(Step2=1;Size;1)
# 606  "spfactor.c"
            for ( Step2 = 1; Step2<=Size; Step2++ ) {
 
                DoRealDirect[Step2] = 0;
            }
# 608  "spfactor.c"
        }
# 606  "spfactor.c"
        Step = 1 + ((Size)>(0) ? (Size) : (0));
# 611  "spfactor.c"
        DoCmplxDirect[Step] = 0;
 
        return ;
# 613  "spfactor.c"
        }
 
    if (!(Mode == 3)) {
# 615  "spfactor.c"
        fflush( (&_iob[1]) );
# 615  "spfactor.c"
        fprintf( (&_iob[2]), "sparse: panic in file `%s' at line %d.\n", "spfactor.c", 615 );
# 615  "spfactor.c"
        fflush( (&_iob[2]) );
# 615  "spfactor.c"
        abort(  );
# 615  "spfactor.c"
        }
# 615  "spfactor.c"
# 595  "spfactor.c"
    ;
# 618  "spfactor.c"
    Nc = (int  *)((*Matrix).MarkowitzRow);
    No = (int  *)((*Matrix).MarkowitzCol);
    Nm = (int  *)((*Matrix).MarkowitzProd);
# 623  "spfactor.c"
    for ( Step = 1; Step<=Size; Step++ ) {
        Nm[Step] = 0;
# 624  "spfactor.c"
        No[Step] = 0;
# 624  "spfactor.c"
        Nc[Step] = No[Step];
 
        pElement = (*Matrix).FirstInCol[Step];
        while ( pElement != 0 ) {
            Nc[Step]++;
            pElement = (*pElement).NextInCol;
# 629  "spfactor.c"
            }
# 632  "spfactor.c"
        pColumn = (*Matrix).FirstInCol[Step];
        while ( (*pColumn).Row < Step ) {
            pElement = (*Matrix).Diag[(*pColumn).Row];
            Nm[Step]++;
            pElement = (*pElement).NextInCol;
# 636  "spfactor.c"
            while ( pElement != 0 ) {
                No[Step]++;
# 636  "spfactor.c"
                pElement = (*pElement).NextInCol;
# 636  "spfactor.c"
                }
 
            pColumn = (*pColumn).NextInCol;
# 638  "spfactor.c"
            }
    }
# 642  "spfactor.c"
    for ( Step = 1; Step<=Size; Step++ ) {
# 680  "spfactor.c"
        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3 * Nc[Step] - 2 * Nm[Step]);
# 683  "spfactor.c"
        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7 * Nc[Step] - 4 * Nm[Step]);
    }
# 696  "spfactor.c"
    return ;
}
# 1  "spfactor.c"
# 723  "spfactor.c"
static int CreateInternalVectors( Matrix )
MatrixPtr Matrix;
 
{
    int Size;
 
# 731  "spfactor.c"
    Size = (*Matrix).Size;
 
    if ((*Matrix).MarkowitzRow == 0) {
        (*Matrix).MarkowitzRow = ((int  *)(malloc( (unsigned int )((4 * (Size + 1))) )));
# 734  "spfactor.c"
        if (((*Matrix).MarkowitzRow) == 0) {
            (*Matrix).Error = 4;
# 735  "spfactor.c"
            }
# 735  "spfactor.c"
        }
 
    if ((*Matrix).MarkowitzCol == 0) {
        (*Matrix).MarkowitzCol = ((int  *)(malloc( (unsigned int )((4 * (Size + 1))) )));
# 738  "spfactor.c"
        if (((*Matrix).MarkowitzCol) == 0) {
            (*Matrix).Error = 4;
# 739  "spfactor.c"
            }
# 739  "spfactor.c"
        }
 
    if ((*Matrix).MarkowitzProd == 0) {
        (*Matrix).MarkowitzProd = ((long  *)(malloc( (unsigned int )((4 * (Size + 2))) )));
# 742  "spfactor.c"
        if (((*Matrix).MarkowitzProd) == 0) {
            (*Matrix).Error = 4;
# 743  "spfactor.c"
            }
# 743  "spfactor.c"
        }
# 748  "spfactor.c"
    if ((*Matrix).DoRealDirect == 0) {
        (*Matrix).DoRealDirect = ((int  *)(malloc( (unsigned int )((4 * (Size + 1))) )));
# 749  "spfactor.c"
        if (((*Matrix).DoRealDirect) == 0) {
            (*Matrix).Error = 4;
# 750  "spfactor.c"
            }
# 750  "spfactor.c"
        }
# 754  "spfactor.c"
    if ((*Matrix).DoCmplxDirect == 0) {
        (*Matrix).DoCmplxDirect = ((int  *)(malloc( (unsigned int )((4 * (Size + 1))) )));
# 755  "spfactor.c"
        if (((*Matrix).DoCmplxDirect) == 0) {
            (*Matrix).Error = 4;
# 756  "spfactor.c"
            }
# 756  "spfactor.c"
        }
# 762  "spfactor.c"
    if ((*Matrix).Intermediate == 0) {
        (*Matrix).Intermediate = ((double  *)(malloc( (unsigned int )((8 * (2 * (Size + 1)))) )));
# 763  "spfactor.c"
        if (((*Matrix).Intermediate) == 0) {
            (*Matrix).Error = 4;
# 764  "spfactor.c"
            }
# 764  "spfactor.c"
        }
# 773  "spfactor.c"
    if ((*Matrix).Error != 4) {
        (*Matrix).InternalVectorsAllocated = 1;
# 774  "spfactor.c"
        }
    return ;
}
# 1  "spfactor.c"
# 812  "spfactor.c"
static int CountMarkowitz( Matrix, RHS, Step )
MatrixPtr Matrix;
    register RealVector RHS;
    int Step;
 
{
    register int Count;
    register int I;
    register int Size = (*Matrix).Size;
    register ElementPtr pElement;
    int ExtRow;
 
# 837  "spfactor.c"
    for ( I = Step; I<=Size; I++ ) {
# 840  "spfactor.c"
        Count = -1;
        pElement = (*Matrix).FirstInRow[I];
        while ( pElement != 0 && (*pElement).Col < Step ) {
            pElement = (*pElement).NextInRow;
# 843  "spfactor.c"
            }
        while ( pElement != 0 ) {
            Count++;
            pElement = (*pElement).NextInRow;
# 846  "spfactor.c"
            }
# 850  "spfactor.c"
        ExtRow = (*Matrix).IntToExtRowMap[I];
# 856  "spfactor.c"
        if (RHS != 0) {
            if ((*Matrix).Complex) {
                if ((RHS[2*ExtRow] != 0.0) || (RHS[2*ExtRow+1] != 0.0)) {
                    Count++;
# 859  "spfactor.c"
                    }
# 859  "spfactor.c"
            } else {
 
                if (RHS[I] != 0.0) {
# 861  "spfactor.c"
                    Count++;
# 861  "spfactor.c"
                    }
# 861  "spfactor.c"
                }
# 861  "spfactor.c"
            }
# 864  "spfactor.c"
        (*Matrix).MarkowitzRow[I] = Count;
    }
# 868  "spfactor.c"
    for ( I = Step; I<=Size; I++ ) {
# 871  "spfactor.c"
        Count = -1;
        pElement = (*Matrix).FirstInCol[I];
        while ( pElement != 0 && (*pElement).Row < Step ) {
            pElement = (*pElement).NextInCol;
# 874  "spfactor.c"
            }
        while ( pElement != 0 ) {
            Count++;
            pElement = (*pElement).NextInCol;
# 877  "spfactor.c"
            }
 
        (*Matrix).MarkowitzCol[I] = Count;
    }
    return ;
}

.