# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by kundert at athos on Fri Jul  1 01:56:11 1988
#
# This archive contains:
#	Makefile	spAllocate.c	spBuild.c	spConfig.h	
#	spDefs.h	spFactor.c	spMatrix.h	spSolve.c	
#

LANG=""; export LANG

echo x - Makefile
cat >Makefile <<'@EOF'
#  Makefile for Sparse1.3
#
#  Ken Kundert
#  UC Berkeley
#
CFLAGS = -O
LINTFLAGS = -lc -lm
SHELL = /bin/sh
CC = cc

HFILES = spConfig.h spDefs.h spMatrix.h
CFILES = spAllocate.c spBuild.c spFactor.c spOutput.c spSolve.c spUtils.c \
	 spFortran.c
OFILES = spAllocate.o spBuild.o spFactor.o spOutput.o spSolve.o spUtils.o \
	 spFortran.o
LIBRARY = sparse.a
DESTINATION = sparse
TESTC = spTest.c
TESTO = spTest.o

SOURCE = $(HFILES) $(CFILES)

$(DESTINATION)	: $(LIBRARY) $(TESTO)
	$(CC) $(CFLAGS) -o $(DESTINATION) $(TESTO) $(LIBRARY) -lm

$(LIBRARY)	: $(OFILES)
	ar r   $(LIBRARY) $?
	ranlib $(LIBRARY)

lint		:
	@lint $(LINTFLAGS) $(CFILES) $(TESTC) | grep -v "but never used"

clean		:
	rm -f $(OFILES) $(LIBRARY) $(TESTO) $(DESTINATION) core

touch		:
	touch -c  $(OFILES) $(LIBRARY)
	ranlib $(LIBRARY)

tags		: $(SOURCE) $(TESTC)
	ctags -t -w $(SOURCE) $(TESTC)

$(OFILES)	: $(HFILES)
$(TESTO)	: $(HFILES)
@EOF

chmod 644 Makefile

echo x - spAllocate.c
cat >spAllocate.c <<'@EOF'
/*
 *  MATRIX ALLOCATION MODULE
 *
 *  Author:                     Advising professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      UC Berkeley
 *
 *  This file contains the allocation and deallocation routines for the
 *  sparse matrix routines.
 *
 *  >>> User accessible functions contained in this file:
 *  spCreate
 *  spDestroy
 *  spError
 *  spWhereSingular
 *  spGetSize
 *  spSetReal
 *  spSetComplex
 *  spFillinCount
 *  spElementCount
 *
 *  >>> Other functions contained in this file:
 *  spcGetElement
 *  InitializeElementBlocks
 *  spcGetFillin
 *  RecordAllocation
 *  AllocateBlockOfAllocationList
 *  EnlargeMatrix
 *  ExpandTranslationArrays
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 */

#ifndef lint
static char copyright[] =
    "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert";
static char RCSid[] =
    "@(#)$Header: spAllocate.c,v 1.2 88/06/18 11:12:22 kundert Exp $";
#endif



/*
 *  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"






/*
 *  MATRIX ALLOCATION
 *
 *  Allocates and initializes the data structures associated with a matrix.
 *
 *  >>> Returned:
 *  A pointer to the matrix is returned cast into the form of a pointer to
 *  a character.  This pointer is then passed and used by the other matrix
 *  routines to refer to a particular matrix.  If an error occurs, the NULL
 *  pointer is returned.
 *
 *  >>> Arguments:
 *  Size  <input>  (int)
 *      Size of matrix or estimate of size of matrix if matrix is EXPANDABLE.
 *  Complex  <input>  (int)
 *      Type of matrix.  If Complex is 0 then the matrix is real, otherwise
 *      the matrix will be complex.  Note that if the routines are not set up
 *      to handle the type of matrix requested, then a spPANIC error will occur.
 *      Further note that if a matrix will be both real and complex, it must
 *      be specified here as being complex.
 *  pError  <output>  (int *)
 *      Returns error flag, needed because function spError() will not work
 *      correctly if spCreate() returns NULL.
 *
 *  >>> Local variables:
 *  AllocatedSize  (int)
 *      The size of the matrix being allocated.
 *  Matrix  (MatrixPtr)
 *      A pointer to the matrix frame being created.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 *  spPANIC
 *  Error is cleared in this routine.
 */

char *
spCreate( Size, Complex, pError )

int  Size, *pError;
BOOLEAN  Complex;
{
register  unsigned  SizePlusOne;
register  MatrixPtr  Matrix;
register  int  I;
int  AllocatedSize;

/* Begin `spCreate'. */
/* Clear error flag. */
    *pError = spOKAY;

/* Test for valid size. */
    if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE))
    {   *pError = spPANIC;
        return NULL;
    }

/* Test for valid type. */
#if NOT spCOMPLEX
    if (Complex)
    {   *pError = spPANIC;
        return NULL;
    }
#endif
#if NOT REAL
    if (NOT Complex)
    {   *pError = spPANIC;
        return NULL;
    }
#endif

/* Create Matrix. */
    AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
    SizePlusOne = (unsigned)(AllocatedSize + 1);

    if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
    {   *pError = spNO_MEMORY;
        return NULL;
    }

/* Initialize matrix */
    Matrix->ID = SPARSE_ID;
    Matrix->Complex = Complex;
    Matrix->PreviousMatrixWasComplex = Complex;
    Matrix->Factored = NO;
    Matrix->Elements = 0;
    Matrix->Error = *pError;
    Matrix->Fillins = 0;
    Matrix->Reordered = NO;
    Matrix->NeedsOrdering = YES;
    Matrix->NumberOfInterchangesIsOdd = NO;
    Matrix->Partitioned = NO;
    Matrix->RowsLinked = NO;
    Matrix->InternalVectorsAllocated = NO;
    Matrix->SingularCol = 0;
    Matrix->SingularRow = 0;
    Matrix->Size = Size;
    Matrix->AllocatedSize = AllocatedSize;
    Matrix->ExtSize = Size;
    Matrix->AllocatedExtSize = AllocatedSize;
    Matrix->CurrentSize = 0;
    Matrix->ExtToIntColMap = NULL;
    Matrix->ExtToIntRowMap = NULL;
    Matrix->IntToExtColMap = NULL;
    Matrix->IntToExtRowMap = NULL;
    Matrix->MarkowitzRow = NULL;
    Matrix->MarkowitzCol = NULL;
    Matrix->MarkowitzProd = NULL;
    Matrix->DoCmplxDirect = NULL;
    Matrix->DoRealDirect = NULL;
    Matrix->Intermediate = NULL;
    Matrix->RelThreshold = DEFAULT_THRESHOLD;
    Matrix->AbsThreshold = 0.0;

    Matrix->TopOfAllocationList = NULL;
    Matrix->RecordsRemaining = 0;
    Matrix->ElementsRemaining = 0;
    Matrix->FillinsRemaining = 0;

    RecordAllocation( Matrix, (char *)Matrix );
    if (Matrix->Error == spNO_MEMORY) goto MemoryError;

/* Take out the trash. */
    Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
    Matrix->TrashCan.Imag = 0.0;
#endif
    Matrix->TrashCan.Row = 0;
    Matrix->TrashCan.Col = 0;
    Matrix->TrashCan.NextInRow = NULL;
    Matrix->TrashCan.NextInCol = NULL;
#if INITIALIZE
    Matrix->TrashCan.pInitInfo = NULL;
#endif

/* Allocate space in memory for Diag pointer vector. */
    CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
    if (Matrix->Diag == NULL)
        goto MemoryError;

/* Allocate space in memory for FirstInCol pointer vector. */
    CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
    if (Matrix->FirstInCol == NULL)
        goto MemoryError;

/* Allocate space in memory for FirstInRow pointer vector. */
    CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
    if (Matrix->FirstInRow == NULL)
        goto MemoryError;

/* Allocate space in memory for IntToExtColMap vector. */
    if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Allocate space in memory for IntToExtRowMap vector. */
    if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Initialize MapIntToExt vectors. */
    for (I = 1; I <= AllocatedSize; I++)
    {   Matrix->IntToExtRowMap[I] = I;
        Matrix->IntToExtColMap[I] = I;
    }

#if TRANSLATE
/* Allocate space in memory for ExtToIntColMap vector. */
    if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Allocate space in memory for ExtToIntRowMap vector. */
    if (( Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Initialize MapExtToInt vectors. */
    for (I = 1; I <= AllocatedSize; I++)
    {   Matrix->ExtToIntColMap[I] = -1;
        Matrix->ExtToIntRowMap[I] = -1;
    }
    Matrix->ExtToIntColMap[0] = 0;
    Matrix->ExtToIntRowMap[0] = 0;
#endif

/* Allocate space for fill-ins and initial set of elements. */
    InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
                                     SPACE_FOR_FILL_INS*AllocatedSize );
    if (Matrix->Error == spNO_MEMORY)
        goto MemoryError;

    return (char *)Matrix;

MemoryError:

/* Deallocate matrix and return no pointer to matrix if there is not enough
   memory. */
    *pError = spNO_MEMORY;
    spDestroy( (char *)Matrix);
    return NULL;
}









/*
 *  ELEMENT ALLOCATION
 *
 *  This routine allocates space for matrix elements. It requests large blocks
 *  of storage from the system and doles out individual elements as required.
 *  This technique, as opposed to allocating elements individually, tends to
 *  speed the allocation process.
 *
 *  >>> Returned:
 *  A pointer to an element.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      A pointer to the first element in the group of elements being allocated.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

ElementPtr
spcGetElement( Matrix )

MatrixPtr Matrix;
{
ElementPtr  pElement;

/* Begin `spcGetElement'. */

/* Allocate block of MatrixElements if necessary. */
    if (Matrix->ElementsRemaining == 0)
    {   pElement = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
        RecordAllocation( Matrix, (char *)pElement );
        if (Matrix->Error == spNO_MEMORY) return NULL;
        Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
        Matrix->NextAvailElement = pElement;
    }

/* Update Element counter and return pointer to Element. */
    Matrix->ElementsRemaining--;
    return Matrix->NextAvailElement++;

}








/*
 *  ELEMENT ALLOCATION INITIALIZATION
 *
 *  This routine allocates space for matrix fill-ins and an initial set of
 *  elements.  Besides being faster than allocating space for elements one
 *  at a time, it tends to keep the fill-ins physically close to the other
 *  matrix elements in the computer memory.  This keeps virtual memory paging
 *  to a minimum.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  InitialNumberOfElements  <input> (int)
 *      This number is used as the size of the block of memory, in
 *      MatrixElements, reserved for elements. If more than this number of
 *      elements are generated, then more space is allocated later.
 *  NumberOfFillinsExpected  <input> (int)
 *      This number is used as the size of the block of memory, in
 *      MatrixElements, reserved for fill-ins. If more than this number of
 *      fill-ins are generated, then more space is allocated, but they may
 *      not be physically close in computer's memory.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      A pointer to the first element in the group of elements being allocated.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static
InitializeElementBlocks( Matrix, InitialNumberOfElements,
                         NumberOfFillinsExpected )

MatrixPtr Matrix;
int  InitialNumberOfElements, NumberOfFillinsExpected;
{
ElementPtr  pElement;

/* Begin `InitializeElementBlocks'. */

/* Allocate block of MatrixElements for elements. */
    pElement = ALLOC(struct MatrixElement, InitialNumberOfElements);
    RecordAllocation( Matrix, (char *)pElement );
    if (Matrix->Error == spNO_MEMORY) return;
    Matrix->ElementsRemaining = InitialNumberOfElements;
    Matrix->NextAvailElement = pElement;

/* Allocate block of MatrixElements for fill-ins. */
    pElement = ALLOC(struct MatrixElement, NumberOfFillinsExpected);
    RecordAllocation( Matrix, (char *)pElement );
    if (Matrix->Error == spNO_MEMORY) return;
    Matrix->FillinsRemaining = NumberOfFillinsExpected;
    Matrix->NextAvailFillin = pElement;

/* Allocate a fill-in list structure. */
    Matrix->FirstFillinListNode = ALLOC(struct FillinListNodeStruct,1);
    RecordAllocation( Matrix, (char *)Matrix->FirstFillinListNode );
    if (Matrix->Error == spNO_MEMORY) return;
    Matrix->LastFillinListNode = Matrix->FirstFillinListNode;

    Matrix->FirstFillinListNode->pFillinList = pElement;
    Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
    Matrix->FirstFillinListNode->Next = NULL;

    return;
}










/*
 *  FILL-IN ALLOCATION
 *
 *  This routine allocates space for matrix fill-ins. It requests large blocks
 *  of storage from the system and doles out individual elements as required.
 *  This technique, as opposed to allocating elements individually, tends to
 *  speed the allocation process.
 *
 *  >>> Returned:
 *  A pointer to the fill-in.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

ElementPtr
spcGetFillin( Matrix )

MatrixPtr Matrix;
{
struct FillinListNodeStruct *pListNode;
ElementPtr  pFillins;

/* Begin `spcGetFillin'. */

#if NOT STRIP OR LINT
    if (Matrix->FillinsRemaining == 0)
        return spcGetElement( Matrix );
#endif
#if STRIP OR LINT

    if (Matrix->FillinsRemaining == 0)
    {   pListNode = Matrix->LastFillinListNode;

/* First see if there are any stripped fill-ins left. */
        if (pListNode->Next != NULL)
        {   Matrix->LastFillinListNode = pListNode = pListNode->Next;
            Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
            Matrix->NextAvailFillin = pListNode->pFillinList;
        }
        else
        {
/* Allocate block of fill-ins. */
            pFillins = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
            RecordAllocation( Matrix, (char *)pFillins );
            if (Matrix->Error == spNO_MEMORY) return NULL;
            Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
            Matrix->NextAvailFillin = pFillins;

/* Allocate a fill-in list structure. */
            pListNode->Next = ALLOC(struct FillinListNodeStruct,1);
            RecordAllocation( Matrix, (char *)pListNode->Next );
            if (Matrix->Error == spNO_MEMORY) return NULL;
            Matrix->LastFillinListNode = pListNode = pListNode->Next;

            pListNode->pFillinList = pFillins;
            pListNode->NumberOfFillinsInList = ELEMENTS_PER_ALLOCATION;
            pListNode->Next = NULL;
        }
    }
#endif

/* Update Fill-in counter and return pointer to Fill-in. */
    Matrix->FillinsRemaining--;
    return Matrix->NextAvailFillin++;
}









/*
 *  RECORD A MEMORY ALLOCATION
 *
 *  This routine is used to record all memory allocations so that the memory
 *  can be freed later.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  AllocatedPtr  <input>  (char *)
 *      The pointer returned by malloc or calloc.  These pointers are saved in
 *      a list so that they can be easily freed.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static
RecordAllocation( Matrix, AllocatedPtr )

MatrixPtr Matrix;
char  *AllocatedPtr;
{
/* Begin `RecordAllocation'. */
/*
 * If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
 * which indicates a spNO_MEMORY error.
 */
    if (AllocatedPtr == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }

/* Allocate block of MatrixElements if necessary. */
    if (Matrix->RecordsRemaining == 0)
    {   AllocateBlockOfAllocationList( Matrix );
        if (Matrix->Error == spNO_MEMORY)
        {   FREE(AllocatedPtr);
            return;
        }
    }

/* Add Allocated pointer to Allocation List. */
    (++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
    Matrix->RecordsRemaining--;
    return;

}








/*
 *  ADD A BLOCK OF SLOTS TO ALLOCATION LIST     
 *
 *  This routine increases the size of the allocation list.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *
 *  >>> Local variables:
 *  ListPtr  (AllocationListPtr)
 *      Pointer to the list that contains the pointers to segments of memory
 *      that were allocated by the operating system for the current matrix.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static
AllocateBlockOfAllocationList( Matrix )

MatrixPtr Matrix;
{
register  int  I;
register  AllocationListPtr  ListPtr;

/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
    ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
    if (ListPtr == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }

/* String entries of allocation list into singly linked list.  List is linked
   such that any record points to the one before it. */

    ListPtr->NextRecord = Matrix->TopOfAllocationList;
    Matrix->TopOfAllocationList = ListPtr;
    ListPtr += ELEMENTS_PER_ALLOCATION;
    for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
    {    ListPtr->NextRecord = ListPtr - 1;
         ListPtr--;
    }

/* Record allocation of space for allocation list on allocation list. */
    Matrix->TopOfAllocationList->AllocatedPtr = (char *)ListPtr;
    Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;

    return;
}








/*
 *  MATRIX DEALLOCATION
 *
 *  Deallocates pointers and elements of Matrix.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to the matrix frame which is to be removed from memory.
 *
 *  >>> Local variables:
 *  ListPtr  (AllocationListPtr)
 *      Pointer into the linked list of pointers to allocated data structures.
 *      Points to pointer to structure to be freed.
 *  NextListPtr  (AllocationListPtr)
 *      Pointer into the linked list of pointers to allocated data structures.
 *      Points to the next pointer to structure to be freed.  This is needed
 *      because the data structure to be freed could include the current node
 *      in the allocation list.
 */

void
spDestroy( eMatrix )

register char *eMatrix;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register  AllocationListPtr  ListPtr, NextListPtr;


/* Begin `spDestroy'. */
    ASSERT( IS_SPARSE( Matrix ) );

/* Deallocate the vectors that are located in the matrix frame. */
    FREE( Matrix->IntToExtColMap );
    FREE( Matrix->IntToExtRowMap );
    FREE( Matrix->ExtToIntColMap );
    FREE( Matrix->ExtToIntRowMap );
    FREE( Matrix->Diag );
    FREE( Matrix->FirstInRow );
    FREE( Matrix->FirstInCol );
    FREE( Matrix->MarkowitzRow );
    FREE( Matrix->MarkowitzCol );
    FREE( Matrix->MarkowitzProd );
    FREE( Matrix->DoCmplxDirect );
    FREE( Matrix->DoRealDirect );
    FREE( Matrix->Intermediate );

/* Sequentially step through the list of allocated pointers freeing pointers
 * along the way. */
    ListPtr = Matrix->TopOfAllocationList;
    while (ListPtr != NULL)
    {   NextListPtr = ListPtr->NextRecord;
        FREE( ListPtr->AllocatedPtr );
        ListPtr = NextListPtr;
    }
    return;
}







/*
 *  RETURN MATRIX ERROR STATUS
 *
 *  This function is used to determine the error status of the given matrix.
 *
 *  >>> Returned:
 *      The error status of the given matrix.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      The matrix for which the error status is desired.
 */

int
spError( eMatrix )

char  *eMatrix;
{
/* Begin `spError'. */

    if (eMatrix != NULL)
    {   ASSERT(((MatrixPtr)eMatrix)->ID == SPARSE_ID);
        return ((MatrixPtr)eMatrix)->Error;
    }
    else return spNO_MEMORY;   /* This error may actually be spPANIC,
                                * no way to tell. */
}









/*
 *  WHERE IS MATRIX SINGULAR
 *
 *  This function returns the row and column number where the matrix was
 *  detected as singular or where a zero was detected on the diagonal.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      The matrix for which the error status is desired.
 *  pRow  <output>  (int *)
 *      The row number.
 *  pCol  <output>  (int *)
 *      The column number.
 */

void
spWhereSingular( eMatrix, pRow, pCol )

char *eMatrix;
int *pRow, *pCol;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;

/* Begin `spWhereSingular'. */
    ASSERT( IS_SPARSE( Matrix ) );

    if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
    {   *pRow = Matrix->SingularRow;
        *pCol = Matrix->SingularCol;
    }
    else *pRow = *pCol = 0;
    return;
}






/*
 *  MATRIX SIZE
 *
 *  Returns the size of the matrix.  Either the internal or external size of
 *  the matrix is returned.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      Pointer to matrix.
 *  External  <input>  (BOOLEAN)
 *      If External is set true, the external size , i.e., the value of the
 *      largest external row or column number encountered is returned.
 *      Otherwise the true size of the matrix is returned.  These two sizes
 *      may differ if the TRANSLATE option is set true.
 */

int
spGetSize( eMatrix, External )

char  *eMatrix;
BOOLEAN  External;
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;

/* Begin `spGetSize'. */
    ASSERT( IS_SPARSE( Matrix ) );

#if TRANSLATE
    if (External)
        return Matrix->ExtSize;
    else
        return Matrix->Size;
#else
    return Matrix->Size;
#endif
}








/*
 *  SET MATRIX COMPLEX OR REAL
 *
 *  Forces matrix to be either real or complex.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      Pointer to matrix.
 */

void
spSetReal( eMatrix )

char *eMatrix;
{
/* Begin `spSetReal'. */

    ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND REAL);
    ((MatrixPtr)eMatrix)->Complex = NO;
    return;
}


void
spSetComplex( eMatrix )

char  *eMatrix;
{
/* Begin `spSetComplex'. */

    ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND spCOMPLEX);
    ((MatrixPtr)eMatrix)->Complex = YES;
    return;
}









/*
 *  ELEMENT OR FILL-IN COUNT
 *
 *  Two functions used to return simple statistics.  Either the number
 *  of total elements, or the number of fill-ins can be returned.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      Pointer to matrix.
 */

int
spFillinCount( eMatrix )

char *eMatrix;
{
/* Begin `spFillinCount'. */

    ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
    return ((MatrixPtr)eMatrix)->Fillins;
}


int
spElementCount( eMatrix )

char  *eMatrix;
{
/* Begin `spElementCount'. */

    ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
    return ((MatrixPtr)eMatrix)->Elements;
}
@EOF

chmod 644 spAllocate.c

echo x - spBuild.c
cat >spBuild.c <<'@EOF'
/*
 *  MATRIX BUILD MODULE
 *
 *  Author:                     Advising professor:
 *     Kenneth S. Kundert           Alberto Sangiovanni-Vincentelli
 *     UC Berkeley
 *
 *  This file contains the routines associated with clearing, loading and
 *  preprocessing the matrix for the sparse matrix routines.
 *
 *  >>> User accessible functions contained in this file:
 *  spClear
 *  spGetElement
 *  spGetAdmittance
 *  spGetQuad
 *  spGetOnes
 *  spInstallInitInfo
 *  spGetInitInfo
 *  spInitialize
 *
 *  >>> Other functions contained in this file:
 *  spcFindElementInCol
 *  Translate
 *  spcCreateElement
 *  spcLinkRows
 *  EnlargeMatrix
 *  ExpandTranslationArrays
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 */

#ifndef lint
static char copyright[] =
    "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert";
static char RCSid[] =
    "@(#)$Header: spBuild.c,v 1.2 88/06/18 11:12:37 kundert Exp $";
#endif




/*
 *  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"






/*
 *  CLEAR MATRIX
 *
 *  Sets every element of the matrix to zero and clears the error flag.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *     Pointer to matrix that is to be cleared.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *     A pointer to the element being cleared.
 */

void
spClear( eMatrix )

char *eMatrix;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement;
register  int  I;

/* Begin `spClear'. */
    ASSERT( IS_SPARSE( Matrix ) );

/* Clear matrix. */
#if spCOMPLEX
    if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex)
    {   for (I = Matrix->Size; I > 0; I--)
        {   pElement = Matrix->FirstInCol[I];
            while (pElement != NULL)
            {   pElement->Real = 0.0;
                pElement->Imag = 0.0;
                pElement = pElement->NextInCol;
            }
        }
    }
    else
#endif
    {   for (I = Matrix->Size; I > 0; I--)
        {   pElement = Matrix->FirstInCol[I];
            while (pElement != NULL)
            {   pElement->Real = 0.0;
                pElement = pElement->NextInCol;
            }
        }
    }

/* Empty the trash. */
    Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
    Matrix->TrashCan.Imag = 0.0;
#endif

    Matrix->Error = spOKAY;
    Matrix->Factored = NO;
    Matrix->SingularCol = 0;
    Matrix->SingularRow = 0;
    Matrix->PreviousMatrixWasComplex = Matrix->Complex;
    return;
}











/*
 *  SINGLE ELEMENT ADDITION TO MATRIX BY INDEX
 *
 *  Finds element [Row,Col] and returns a pointer to it.  If element is
 *  not found then it is created and spliced into matrix.  This routine
 *  is only to be used after spCreate() and before spMNA_Preorder(),
 *  spFactor() or spOrderAndFactor().  Returns a pointer to the
 *  Real portion of a MatrixElement.  This pointer is later used by
 *  spADD_xxx_ELEMENT to directly access element.
 *
 *  >>> Returns:
 *  Returns a pointer to the element.  This pointer is then used to directly
 *  access the element during successive builds.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *     Pointer to the matrix that the element is to be added to.
 *  Row  <input>  (int)
 *     Row index for element.  Must be in the range of [0..Size] unless
 *     the options EXPANDABLE or TRANSLATE are used. Elements placed in
 *     row zero are discarded.  In no case may Row be less than zero.
 *  Col  <input>  (int)
 *     Column index for element.  Must be in the range of [0..Size] unless
 *     the options EXPANDABLE or TRANSLATE are used. Elements placed in
 *     column zero are discarded.  In no case may Col be less than zero.
 *
 *  >>> Local variables:
 *  pElement  (RealNumber *)
 *     Pointer to the element.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 *  Error is not cleared in this routine.
 */

RealNumber *
spGetElement( eMatrix, Row, Col )

char *eMatrix;
int  Row, Col;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
RealNumber  *pElement;
ElementPtr spcFindElementInCol();
void  Translate();

/* Begin `spGetElement'. */
    ASSERT( IS_SPARSE( Matrix ) AND Row >= 0 AND Col >= 0 );

    if ((Row == 0) OR (Col == 0))
        return &Matrix->TrashCan.Real;

#if NOT TRANSLATE
    ASSERT(Matrix->NeedsOrdering);
#endif

#if TRANSLATE
    Translate( Matrix, &Row, &Col );
    if (Matrix->Error == spNO_MEMORY) return NULL;
#endif

#if NOT TRANSLATE
#if NOT EXPANDABLE
    ASSERT(Row <= Matrix->Size AND Col <= Matrix->Size);
#endif

#if EXPANDABLE
/* Re-size Matrix if necessary. */
    if ((Row > Matrix->Size) OR (Col > Matrix->Size))
        EnlargeMatrix( Matrix, MAX(Row, Col) );
    if (Matrix->Error == spNO_MEMORY) return NULL;
#endif
#endif

/*
 * The condition part of the following if statement tests to see if the
 * element resides along the diagonal, if it does then it tests to see
 * if the element has been created yet (Diag pointer not NULL).  The
 * pointer to the element is then assigned to Element after it is cast
 * into a pointer to a RealNumber.  This casting makes the pointer into
 * a pointer to Real.  This statement depends on the fact that Real
 * is the first record in the MatrixElement structure.
 */

    if ((Row != Col) OR ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL))
    {
/*
 * Element does not exist or does not reside along diagonal.  Search
 * column for element.  As in the if statement above, the pointer to the
 * element which is returned by spcFindElementInCol is cast into a
 * pointer to Real, a RealNumber.
 */
        pElement = (RealNumber*)spcFindElementInCol( Matrix,
                                                     &(Matrix->FirstInCol[Col]),
                                                     Row, Col, YES );
    }
    return pElement;
}











/*
 *  FIND ELEMENT BY SEARCHING COLUMN
 *
 *  Searches column starting at element specified at PtrAddr and finds element
 *  in Row. If Element does not exists, it is created. The pointer to the
 *  element is returned.
 *
 *  >>> Returned:
 *  A pointer to the desired element:
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to Matrix.
 *  LastAddr  <input-output>  (ElementPtr *)
 *      Address of pointer that initially points to the element in Col at which
 *      the search is started.  The pointer in this location may be changed if
 *      a fill-in is required in and adjacent element. For this reason it is
 *      important that LastAddr be the address of a FirstInCol or a NextInCol
 *      rather than a temporary variable.
 *  Row  <input>  (int)
 *      Row being searched for.
 *  Col  (int)
 *      Column being searched.
 *  CreateIfMissing  <input>  (BOOLEAN)
 *      Indicates what to do if element is not found, create one or return a
 *      NULL pointer.
 *
 *  Local variables:
 *  pElement  (ElementPtr)
 *      Pointer used to search through matrix.
 */

ElementPtr
spcFindElementInCol( Matrix, LastAddr, Row, Col, CreateIfMissing )

MatrixPtr Matrix;
register ElementPtr *LastAddr;
register int  Row;
int  Col;
BOOLEAN  CreateIfMissing;
{
register  ElementPtr  pElement;
ElementPtr  spcCreateElement();

/* Begin `spcFindElementInCol'. */
    pElement = *LastAddr;

/* Search for element. */
    while (pElement != NULL)
    {   if (pElement->Row < Row)
        {
/* Have not reached element yet. */
            LastAddr = &(pElement->NextInCol);
            pElement = pElement->NextInCol;
        }
        else if (pElement->Row == Row)
        {
/* Reached element. */
            return pElement;
        }
        else break;  /* while loop */
    }

/* Element does not exist and must be created. */
    if (CreateIfMissing)
       return spcCreateElement( Matrix, Row, Col, LastAddr, NO );
    else return NULL;
}








#if TRANSLATE

/*
 *  TRANSLATE EXTERNAL INDICES TO INTERNAL
 *
 *  Convert internal row and column numbers to internal row and column numbers.
 *  Also updates Ext/Int maps.
 *
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  Row  <input/output>  (int *)
 *     Upon entry Row is either a external row number of an external node
 *     number.  Upon entry, the internal equivalent is supplied.
 *  Col  <input/output>  (int *)
 *     Upon entry Column is either a external column number of an external node
 *     number.  Upon entry, the internal equivalent is supplied.
 *
 *  >>> Local variables:
 *  ExtCol  (int)
 *     Temporary variable used to hold the external column or node number
 *     during the external to internal column number translation.
 *  ExtRow  (int)
 *     Temporary variable used to hold the external row or node number during
 *     the external to internal row number translation.
 *  IntCol  (int)
 *     Temporary variable used to hold the internal column or node number
 *     during the external to internal column number translation.
 *  IntRow  (int)
 *     Temporary variable used to hold the internal row or node number during
 *     the external to internal row number translation.
 */

static void
Translate( Matrix, Row, Col )

MatrixPtr Matrix;
int  *Row, *Col;
{
register int IntRow, IntCol, ExtRow, ExtCol;

/* Begin `Translate'. */
    ExtRow = *Row;
    ExtCol = *Col;

/* Expand translation arrays if necessary. */
    if ((ExtRow > Matrix->AllocatedExtSize) OR
        (ExtCol > Matrix->AllocatedExtSize))
    {
        ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) );
        if (Matrix->Error == spNO_MEMORY) return;
    }

/* Set ExtSize if necessary. */
    if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize))
        Matrix->ExtSize = MAX(ExtRow, ExtCol);

/* Translate external row or node number to internal row or node number. */
    if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1)
    {   Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
        Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize;
        IntRow = Matrix->CurrentSize;

#if NOT EXPANDABLE
        ASSERT(IntRow <= Matrix->Size);
#endif

#if EXPANDABLE
/* Re-size Matrix if necessary. */
        if (IntRow > Matrix->Size)
            EnlargeMatrix( Matrix, IntRow );
        if (Matrix->Error == spNO_MEMORY) return;
#endif

        Matrix->IntToExtRowMap[IntRow] = ExtRow;
        Matrix->IntToExtColMap[IntRow] = ExtRow;
    }

/* Translate external column or node number to internal column or node number.*/
    if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1)
    {   Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize;
        Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize;
        IntCol = Matrix->CurrentSize;

#if NOT EXPANDABLE
        ASSERT(IntCol <= Matrix->Size);
#endif

#if EXPANDABLE
/* Re-size Matrix if necessary. */
        if (IntCol > Matrix->Size)
            EnlargeMatrix( Matrix, IntCol );
        if (Matrix->Error == spNO_MEMORY) return;
#endif

        Matrix->IntToExtRowMap[IntCol] = ExtCol;
        Matrix->IntToExtColMap[IntCol] = ExtCol;
    }

    *Row = IntRow;
    *Col = IntCol;
    return;
}
#endif






#if QUAD_ELEMENT
/*
 *  ADDITION OF ADMITTANCE TO MATRIX BY INDEX
 *
 *  Performs same function as spGetElement except rather than one
 *  element, all four Matrix elements for a floating component are
 *  added.  This routine also works if component is grounded.  Positive
 *  elements are placed at [Node1,Node2] and [Node2,Node1].  This
 *  routine is only to be used after spCreate() and before
 *  spMNA_Preorder(), spFactor() or spOrderAndFactor().
 *
 *  >>> Returns:
 *  Error code.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *     Pointer to the matrix that component is to be entered in.
 *  Node1  <input>  (int)
 *     Row and column indices for elements. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Node zero is the
 *     ground node.  In no case may Node1 be less than zero.
 *  Node2  <input>  (int)
 *     Row and column indices for elements. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Node zero is the
 *     ground node.  In no case may Node2 be less than zero.
 *  Template  <output>  (struct spTemplate *)
 *     Collection of pointers to four elements that are later used to directly
 *     address elements.  User must supply the template, this routine will
 *     fill it.
 *
 *  Possible errors:
 *  spNO_MEMORY
 *  Error is not cleared in this routine.
 */

int
spGetAdmittance( Matrix, Node1, Node2, Template )

char  *Matrix;
int  Node1, Node2;
struct  spTemplate  *Template;
{

/* Begin `spGetAdmittance'. */
    Template->Element1 = spGetElement(Matrix, Node1, Node1 );
    Template->Element2 = spGetElement(Matrix, Node2, Node2 );
    Template->Element3Negated = spGetElement( Matrix, Node2, Node1 );
    Template->Element4Negated = spGetElement( Matrix, Node1, Node2 );
    if
    (   (Template->Element1 == NULL)
        OR (Template->Element2 == NULL)
        OR (Template->Element3Negated == NULL)
        OR (Template->Element4Negated == NULL)
    )   return spNO_MEMORY;

    if (Node1 == 0)
        SWAP( RealNumber*, Template->Element1, Template->Element2 );

    return spOKAY;
}
#endif /* QUAD_ELEMENT */









#if QUAD_ELEMENT
/*
 *  ADDITION OF FOUR ELEMENTS TO MATRIX BY INDEX
 *
 *  Similar to spGetAdmittance, except that spGetAdmittance only
 *  handles 2-terminal components, whereas spGetQuad handles simple
 *  4-terminals as well.  These 4-terminals are simply generalized
 *  2-terminals with the option of having the sense terminals different
 *  from the source and sink terminals.  spGetQuad adds four
 *  elements to the matrix.  Positive elements occur at Row1,Col1
 *  Row2,Col2 while negative elements occur at Row1,Col2 and Row2,Col1.
 *  The routine works fine if any of the rows and columns are zero.
 *  This routine is only to be used after spCreate() and before
 *  spMNA_Preorder(), spFactor() or spOrderAndFactor()
 *  unless TRANSLATE is set true.
 *
 *  >>> Returns:
 *  Error code.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *     Pointer to the matrix that component is to be entered in.
 *  Row1  <input>  (int)
 *     First row index for elements. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground row.  In no case may Row1 be less than zero.
 *  Row2  <input>  (int)
 *     Second row index for elements. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground row.  In no case may Row2 be less than zero.
 *  Col1  <input>  (int)
 *     First column index for elements. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground column.  In no case may Col1 be less than zero.
 *  Col2  <input>  (int)
 *     Second column index for elements. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground column.  In no case may Col2 be less than zero.
 *  Template  <output>  (struct spTemplate *)
 *     Collection of pointers to four elements that are later used to directly
 *     address elements.  User must supply the template, this routine will
 *     fill it.
 *  Real  <input>  (RealNumber)
 *     Real data to be added to elements.
 *  Imag  <input>  (RealNumber)
 *     Imag data to be added to elements.  If matrix is real, this argument
 *     may be deleted.
 *
 *  Possible errors:
 *  spNO_MEMORY
 *  Error is not cleared in this routine.
 */

int
spGetQuad( Matrix, Row1, Row2, Col1, Col2, Template )

char  *Matrix;
int  Row1, Row2, Col1, Col2;
struct  spTemplate  *Template;
{
/* Begin `spGetQuad'. */
    Template->Element1 = spGetElement( Matrix, Row1, Col1);
    Template->Element2 = spGetElement( Matrix, Row2, Col2 );
    Template->Element3Negated = spGetElement( Matrix, Row2, Col1 );
    Template->Element4Negated = spGetElement( Matrix, Row1, Col2 );
    if
    (   (Template->Element1 == NULL)
        OR (Template->Element2 == NULL)
        OR (Template->Element3Negated == NULL)
        OR (Template->Element4Negated == NULL)
    )   return spNO_MEMORY;

    if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real)
        SWAP( RealNumber *, Template->Element1, Template->Element2 );

    return spOKAY;
}
#endif /* QUAD_ELEMENT */









#if QUAD_ELEMENT
/*
 *  ADDITION OF FOUR STRUCTURAL ONES TO MATRIX BY INDEX
 *
 *  Performs similar function to spGetQuad() except this routine is
 *  meant for components that do not have an admittance representation.
 *
 *  The following stamp is used:
 *         Pos  Neg  Eqn
 *  Pos  [  .    .    1  ]
 *  Neg  [  .    .   -1  ]
 *  Eqn  [  1   -1    .  ]
 *
 *  >>> Returns:
 *  Error code.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *     Pointer to the matrix that component is to be entered in.
 *  Pos  <input>  (int)
 *     See stamp above. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground row.  In no case may Pos be less than zero.
 *  Neg  <input>  (int)
 *     See stamp above. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground row.  In no case may Neg be less than zero.
 *  Eqn  <input>  (int)
 *     See stamp above. Must be in the range of [0..Size]
 *     unless the options EXPANDABLE or TRANSLATE are used. Zero is the
 *     ground row.  In no case may Eqn be less than zero.
 *  Template  <output>  (struct spTemplate *)
 *     Collection of pointers to four elements that are later used to directly
 *     address elements.  User must supply the template, this routine will
 *     fill it.
 *
 *  Possible errors:
 *  spNO_MEMORY
 *  Error is not cleared in this routine.
 */

int
spGetOnes(Matrix, Pos, Neg, Eqn, Template)

char  *Matrix;
int  Pos, Neg, Eqn;
struct  spTemplate  *Template;
{
/* Begin `spGetOnes'. */
    Template->Element4Negated = spGetElement( Matrix, Neg, Eqn );
    Template->Element3Negated = spGetElement( Matrix, Eqn, Neg );
    Template->Element2 = spGetElement( Matrix, Pos, Eqn );
    Template->Element1 = spGetElement( Matrix, Eqn, Pos );
    if
    (   (Template->Element1 == NULL)
        OR (Template->Element2 == NULL)
        OR (Template->Element3Negated == NULL)
        OR (Template->Element4Negated == NULL)
    )   return spNO_MEMORY;

    spADD_REAL_QUAD( *Template, 1.0 );
    return spOKAY;
}
#endif /* QUAD_ELEMENT */







/*
 *
 *  CREATE AND SPLICE ELEMENT INTO MATRIX
 *
 *  This routine is used to create new matrix elements and splice them into the
 *  matrix.
 *
 *  >>> Returned:
 *  A pointer to the element that was created is returned.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Row  <input>  (int)
 *      Row index for element.
 *  Col  <input>  (int)
 *      Column index for element.
 *  LastAddr  <input-output>  (ElementPtr *)
 *      This contains the address of the pointer to the element just above the
 *      one being created. It is used to speed the search and it is updated with
 *      address of the created element.
 *  Fillin  <input>  (BOOLEAN)
 *      Flag that indicates if created element is to be a fill-in.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      Pointer to an element in the matrix. It is used to refer to the newly
 *      created element and to restring the pointers of the element's row and
 *      column.
 *  pLastElement  (ElementPtr)
 *      Pointer to the element in the matrix that was just previously pointed
 *      to by pElement. It is used to restring the pointers of the element's
 *      row and column.
 *  pCreatedElement  (ElementPtr)
 *      Pointer to the desired element, the one that was just created.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

ElementPtr
spcCreateElement( Matrix, Row, Col, LastAddr, Fillin )

MatrixPtr Matrix;
int  Row;
register int  Col;
register ElementPtr  *LastAddr;
BOOLEAN Fillin;
{
register  ElementPtr  pElement, pLastElement;
ElementPtr  pCreatedElement, spcGetElement(), spcGetFillin();

/* Begin `spcCreateElement'. */

    if (Matrix->RowsLinked)
    {
/* Row pointers cannot be ignored. */
        if (Fillin)
        {   pElement = spcGetFillin( Matrix );
            Matrix->Fillins++;
        }
        else
        {   pElement = spcGetElement( Matrix );
            Matrix->NeedsOrdering = YES;
        }
        if (pElement == NULL) return NULL;

/* If element is on diagonal, store pointer in Diag. */
        if (Row == Col) Matrix->Diag[Row] = pElement;

/* Initialize Element. */
        pCreatedElement = pElement;
        pElement->Row = Row;
        pElement->Col = Col;
        pElement->Real = 0.0;
#if spCOMPLEX
        pElement->Imag = 0.0;
#endif
#if INITIALIZE
        pElement->pInitInfo = NULL;
#endif

/* Splice element into column. */
        pElement->NextInCol = *LastAddr;
        *LastAddr = pElement;

 /* Search row for proper element position. */
        pElement = Matrix->FirstInRow[Row];
        pLastElement = NULL;
        while (pElement != NULL)
        {
/* Search for element row position. */
            if (pElement->Col < Col)
            {
/* Have not reached desired element. */
                pLastElement = pElement;
                pElement = pElement->NextInRow;
            }
            else pElement = NULL;
        }

/* Splice element into row. */
        pElement = pCreatedElement;
        if (pLastElement == NULL)
        {
/* Element is first in row. */
            pElement->NextInRow = Matrix->FirstInRow[Row];
            Matrix->FirstInRow[Row] = pElement;
        }
        else
/* Element is not first in row. */
        {
            pElement->NextInRow = pLastElement->NextInRow;
            pLastElement->NextInRow = pElement;
        }

    }
    else
    {
/*
 * Matrix has not been factored yet.  Thus get element rather than fill-in.
 * Also, row pointers can be ignored.
 */

/* Allocate memory for Element. */
        pElement = spcGetElement( Matrix );
        if (pElement == NULL) return NULL;

/* If element is on diagonal, store pointer in Diag. */
        if (Row == Col) Matrix->Diag[Row] = pElement;

/* Initialize Element. */
        pCreatedElement = pElement;
        pElement->Row = Row;
#if DEBUG
        pElement->Col = Col;
#endif
        pElement->Real = 0.0;
#if spCOMPLEX
        pElement->Imag = 0.0;
#endif
#if INITIALIZE
        pElement->pInitInfo = NULL;
#endif

/* Splice element into column. */
        pElement->NextInCol = *LastAddr;
        *LastAddr = pElement;
    }

    Matrix->Elements++;
    return pCreatedElement;
}








/*
 *
 *  LINK ROWS
 *
 *  This routine is used to generate the row links.  The spGetElement()
 *  routines do not create row links, which are needed by the spFactor()
 *  routines.
 *
 *  >>>  Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      Pointer to an element in the matrix.
 *  FirstInRowEntry  (ElementPtr *)
 *      A pointer into the FirstInRow array.  Points to the FirstInRow entry
 *      currently being operated upon.
 *  FirstInRowArray  (ArrayOfElementPtrs)
 *      A pointer to the FirstInRow array.  Same as Matrix->FirstInRow but
 *      resides in a register and requires less indirection so is faster to
 *      use.
 *  Col  (int)
 *      Column currently being operated upon.
 */

spcLinkRows( Matrix )

MatrixPtr Matrix;
{
register  ElementPtr  pElement, *FirstInRowEntry;
register  ArrayOfElementPtrs  FirstInRowArray;
register  int  Col;

/* Begin `spcLinkRows'. */
    FirstInRowArray = Matrix->FirstInRow;
    for (Col = Matrix->Size; Col >= 1; Col--)
    {
/* Generate row links for the elements in the Col'th column. */
        pElement = Matrix->FirstInCol[Col];

        while (pElement != NULL)
        {   pElement->Col = Col;
            FirstInRowEntry = &FirstInRowArray[pElement->Row];
            pElement->NextInRow = *FirstInRowEntry;
            *FirstInRowEntry = pElement;
            pElement = pElement->NextInCol;
        }
    }
    Matrix->RowsLinked = YES;
    return;
}








/*
 *  ENLARGE MATRIX
 *
 *  Increases the size of the matrix.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  NewSize  <input>  (int)
 *     The new size of the matrix.
 *
 *  >>> Local variables:
 *  OldAllocatedSize  (int)
 *     The allocated size of the matrix before it is expanded.
 */

static
EnlargeMatrix( Matrix, NewSize )

MatrixPtr Matrix;
register int  NewSize;
{
register int I, OldAllocatedSize = Matrix->AllocatedSize;

/* Begin `EnlargeMatrix'. */
    Matrix->Size = NewSize;

    if (NewSize <= OldAllocatedSize)
        return;

/* Expand the matrix frame. */
    NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
    Matrix->AllocatedSize = NewSize;

    if (( REALLOC(Matrix->IntToExtColMap, int, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }
    if (( REALLOC(Matrix->IntToExtRowMap, int, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }
    if (( REALLOC(Matrix->Diag, ElementPtr, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }
    if (( REALLOC(Matrix->FirstInCol, ElementPtr, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }
    if (( REALLOC(Matrix->FirstInRow, ElementPtr, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }

/*
 * Destroy the Markowitz and Intermediate vectors, they will be recreated
 * in spOrderAndFactor().
 */
    FREE( Matrix->MarkowitzRow );
    FREE( Matrix->MarkowitzCol );
    FREE( Matrix->MarkowitzProd );
    FREE( Matrix->DoRealDirect );
    FREE( Matrix->DoCmplxDirect );
    FREE( Matrix->Intermediate );
    Matrix->InternalVectorsAllocated = NO;

/* Initialize the new portion of the vectors. */
    for (I = OldAllocatedSize+1; I <= NewSize; I++)
    {   Matrix->IntToExtColMap[I] = I;
        Matrix->IntToExtRowMap[I] = I;
        Matrix->Diag[I] = NULL;
        Matrix->FirstInRow[I] = NULL;
        Matrix->FirstInCol[I] = NULL;
    }

    return;
}








#if TRANSLATE

/*
 *  EXPAND TRANSLATION ARRAYS
 *
 *  Increases the size arrays that are used to translate external to internal
 *  row and column numbers.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  NewSize  <input>  (int)
 *     The new size of the translation arrays.
 *
 *  >>> Local variables:
 *  OldAllocatedSize  (int)
 *     The allocated size of the translation arrays before being expanded.
 */

static
ExpandTranslationArrays( Matrix, NewSize )

MatrixPtr Matrix;
register int  NewSize;
{
register int I, OldAllocatedSize = Matrix->AllocatedExtSize;

/* Begin `ExpandTranslationArrays'. */
    Matrix->ExtSize = NewSize;

    if (NewSize <= OldAllocatedSize)
        return;

/* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */
    NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
    Matrix->AllocatedExtSize = NewSize;

    if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }
    if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL)
    {   Matrix->Error = spNO_MEMORY;
        return;
    }

/* Initialize the new portion of the vectors. */
    for (I = OldAllocatedSize+1; I <= NewSize; I++)
    {   Matrix->ExtToIntRowMap[I] = -1;
        Matrix->ExtToIntColMap[I] = -1;
    }

    return;
}
#endif









#if INITIALIZE
/*
 *   INITIALIZE MATRIX
 *
 *   With the INITIALIZE compiler option (see spConfig.h) set true,
 *   Sparse allows the user to keep initialization information with each
 *   structurally nonzero matrix element.  Each element has a pointer
 *   that is set and used by the user.  The user can set this pointer
 *   using spInstallInitInfo and may be read using spGetInitInfo.  Both
 *   may be used only after the element exists.  The function
 *   spInitialize() is a user customizable way to initialize the matrix.
 *   Passed to this routine is a function pointer.  spInitialize() sweeps
 *   through every element in the matrix and checks the pInitInfo
 *   pointer (the user supplied pointer).  If the pInitInfo is NULL,
 *   which is true unless the user changes it (almost always true for
 *   fill-ins), then the element is zeroed.  Otherwise, the function
 *   pointer is called and passed the pInitInfo pointer as well as the
 *   element pointer and the external row and column numbers.  If the
 *   user sets the value of each element, then spInitialize() replaces
 *   spClear().
 *
 *   The user function is expected to return a nonzero integer if there
 *   is a fatal error and zero otherwise.  Upon encountering a nonzero
 *   return code, spInitialize() terminates and returns the error code.
 *
 *   >>> Arguments:
 *   Matrix  <input>  (char *)
 *       Pointer to matrix.
 *
 *   >>> Possible Errors:
 *   Returns nonzero if error, zero otherwise.
 */

void
spInstallInitInfo( pElement, pInitInfo )

RealNumber *pElement;
char *pInitInfo;
{
/* Begin `spInstallInitInfo'. */
    ASSERT(pElement != NULL);

    ((ElementPtr)pElement)->pInitInfo = pInitInfo;
}


char *
spGetInitInfo( pElement )

RealNumber *pElement;
{
/* Begin `spGetInitInfo'. */
    ASSERT(pElement != NULL);

    return (char *)((ElementPtr)pElement)->pInitInfo;
}


int
spInitialize( eMatrix, pInit )

char *eMatrix;
int (*pInit)();
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register ElementPtr pElement;
int J, Error, Col;

/* Begin `spInitialize'. */
    ASSERT( IS_SPARSE( Matrix ) );

#if spCOMPLEX
/* Clear imaginary part of matrix if matrix is real but was complex. */
    if (Matrix->PreviousMatrixWasComplex AND NOT Matrix->Complex)
    {   for (J = Matrix->Size; J > 0; J--)
        {   pElement = Matrix->FirstInCol[J];
            while (pElement != NULL)
            {   pElement->Imag = 0.0;
                pElement = pElement->NextInCol;
            }
        }
    }
#endif /* spCOMPLEX */

/* Initialize the matrix. */
    for (J = Matrix->Size; J > 0; J--)
    {   pElement = Matrix->FirstInCol[J];
        Col = Matrix->IntToExtColMap[J];
        while (pElement != NULL)
        {   if (pElement->pInitInfo == NULL)
            {   pElement->Real = 0.0;
#               if spCOMPLEX
                    pElement->Imag = 0.0;
#               endif
            }
            else
            {   Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo,
                                 Matrix->IntToExtRowMap[pElement->Row], Col);
                if (Error)
                {   Matrix->Error = spFATAL;
                    return Error;
                }

            }
            pElement = pElement->NextInCol;
        }
    }

/* Empty the trash. */
    Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
    Matrix->TrashCan.Imag = 0.0;
#endif

    Matrix->Error = spOKAY;
    Matrix->Factored = NO;
    Matrix->SingularCol = 0;
    Matrix->SingularRow = 0;
    Matrix->PreviousMatrixWasComplex = Matrix->Complex;
    return 0;
}
#endif /* INITIALIZE */
@EOF

chmod 644 spBuild.c

echo x - spConfig.h
cat >spConfig.h <<'@EOF'
/*
 *  CONFIGURATION MACRO DEFINITIONS for sparse matrix routines
 *
 *  Author:                     Advising professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      U.C. Berkeley
 *
 *  This file contains macros for the sparse matrix routines that are used
 *  to define the personality of the routines.  The user is expected to
 *  modify this file to maximize the performance of the routines with
 *  his/her matrices.
 *
 *  Macros are distinguished by using solely capital letters in their
 *  identifiers.  This contrasts with C defined identifiers which are
 *  strictly lower case, and program variable and procedure names which use
 *  both upper and lower case.
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 *
 *  $Date: 88/06/18 11:13:29 $
 *  $Revision: 1.2 $
 */


#ifndef spCONFIG_DEFS
#define spCONFIG_DEFS




#ifdef spINSIDE_SPARSE
/*
 *  OPTIONS
 *
 *  These are compiler options.  Set each option to one to compile that
 *  section of the code.  If a feature is not desired, set the macro
 *  to NO. Recommendations are given in brackets, [ignore them].
 *
 *  >>> Option descriptions:
 *  Arithmetic Precision
 *      The precision of the arithmetic used by Sparse can be set by
 *      changing changing the spREAL macro.  This macro is
 *      contained in the file spMatrix.h.  It is strongly suggested to
 *      used double precision with circuit simulators.  Note that
 *      because C always performs arithmetic operations in double
 *      precision, the only benefit to using single precision is that
 *      less storage is required.  There is often a noticeable speed
 *      penalty when using single precision.  Sparse internally refers
 *      to a spREAL as a RealNumber.
 *  REAL
 *      This specifies that the routines are expected to handle real
 *      systems of equations.  The routines can be compiled to handle
 *      both real and complex systems at the same time, but there is a
 *      slight speed and memory advantage if the routines are complied
 *      to handle only real systems of equations.
 *  spCOMPLEX
 *      This specifies that the routines will be complied to handle
 *      complex systems of equations.
 *  EXPANDABLE
 *      Setting this compiler flag true (1) makes the matrix
 *      expandable before it has been factored.  If the matrix is
 *      expandable, then if an element is added that would be
 *      considered out of bounds in the current matrix, the size of
 *      the matrix is increased to hold that element.  As a result,
 *      the size of the matrix need not be known before the matrix is
 *      built.  The matrix can be allocated with size zero and
 *      expanded.
 *  TRANSLATE
 *      This option allows the set of external row and column numbers
 *      to be non-packed.  In other words, the row and column numbers
 *      do not have to be contiguous.  The priced paid for this
 *      flexibility is that when TRANSLATE is set true, the time
 *      required to initially build the matrix will be greater because
 *      the external row and column number must be translated into
 *      internal equivalents.  This translation brings about other
 *      benefits though.  First, the spGetElement() and
 *      spGetAdmittance() routines may be used after the matrix has
 *      been factored.  Further, elements, and even rows and columns,
 *      may be added to the matrix, and row and columns may be deleted
 *      from the matrix, after it has been factored.  Note that when
 *      the set of row and column number is not a packed set, neither
 *      are the RHS and Solution vectors.  Thus the size of these
 *      vectors must be at least as large as the external size, which
 *      is the value of the largest given row or column numbers.
 *  INITIALIZE
 *      Causes the spInitialize(), spGetInitInfo(), and
 *      spInstallInitInfo() routines to be compiled.  These routines
 *      allow the user to store and read one pointer in each nonzero
 *      element in the matrix.  spInitialize() then calls a user
 *      specified function for each structural nonzero in the matrix,
 *      and includes this pointer as well as the external row and
 *      column numbers as arguments.  This allows the user to write
 *      custom matrix initialization routines.
 *  DIAGONAL_PIVOTING
 *      Many matrices, and in particular node- and modified-node
 *      admittance matrices, tend to be nearly symmetric and nearly
 *      diagonally dominant.  For these matrices, it is a good idea to
 *      select pivots from the diagonal.  With this option enabled,
 *      this is exactly what happens, though if no satisfactory pivot
 *      can be found on the diagonal, an off-diagonal pivot will be
 *      used.  If this option is disabled, Sparse does not
 *      preferentially search the diagonal.  Because of this, Sparse
 *      has a wider variety of pivot candidates available, and so
 *      presumably fewer fill-ins will be created.  However, the
 *      initial pivot selection process will take considerably longer.
 *      If working with node admittance matrices, or other matrices
 *      with a strong diagonal, it is probably best to use
 *      DIAGONAL_PIVOTING for two reasons.  First, accuracy will be
 *      better because pivots will be chosen from the large diagonal
 *      elements, thus reducing the chance of growth.  Second, a near
 *      optimal ordering will be chosen quickly.  If the class of
 *      matrices you are working with does not have a strong diagonal,
 *      do not use DIAGONAL_PIVOTING, but consider using a larger
 *      threshold.  When DIAGONAL_PIVOTING is turned off, the following
 *      options and constants are not used: MODIFIED_MARKOWITZ,
 *      MAX_MARKOWITZ_TIES, and TIES_MULTIPLIER.
 *  ARRAY_OFFSET
 *      This determines whether arrays start at an index of zero or one.
 *      This option is necessitated by the fact that standard C
 *      convention dictates that arrays begin with an index of zero but
 *      the standard mathematic convention states that arrays begin with
 *      an index of one.  So if you prefer to start your arrays with
 *      zero, or your calling Sparse from FORTRAN, set ARRAY_OFFSET to
 *      NO or 0.  Otherwise, set ARRAY_OFFSET to YES or 1.  Note that if
 *      you use an offset of one, the arrays that you pass to Sparse
 *      must have an allocated length of one plus the size of the
 *      matrix.  ARRAY_OFFSET must be either 0 or 1, no other offsets
 *      are valid.
 *  spSEPARATED_COMPLEX_VECTORS
 *      This specifies the format for complex vectors.  If this is set
 *      false then a complex vector is made up of one double sized
 *      array of RealNumber's in which the real and imaginary numbers
 *      are placed in the alternately array in the array.  In other
 *      words, the first entry would be Complex[1].Real, then comes
 *      Complex[1].Imag, then Complex[1].Real, etc.  If
 *      spSEPARATED_COMPLEX_VECTORS is set true, then each complex
 *      vector is represented by two arrays of RealNumbers, one with
 *      the real terms, the other with the imaginary. [NO]
 *  MODIFIED_MARKOWITZ
 *      This specifies that the modified Markowitz method of pivot
 *      selection is to be used.  The modified Markowitz method differs
 *      from standard Markowitz in two ways.  First, under modified
 *      Markowitz, the search for a pivot can be terminated early if a
 *      adequate (in terms of sparsity) pivot candidate is found.
 *      Thus, when using modified Markowitz, the initial factorization
 *      can be faster, but at the expense of a suboptimal pivoting
 *      order that may slow subsequent factorizations.  The second
 *      difference is in the way modified Markowitz breaks Markowitz
 *      ties.  When two or more elements are pivot candidates and they
 *      all have the same Markowitz product, then the tie is broken by
 *      choosing the element that is best numerically.  The numerically
 *      best element is the one with the largest ratio of its magnitude
 *      to the magnitude of the largest element in the same column,
 *      excluding itself.  The modified Markowitz method results in
 *      marginally better accuracy.  This option is most appropriate
 *      for use when working with very large matrices where the initial
 *      factor time represents an unacceptable burden. [NO]
 *  DELETE
 *      This specifies that the spDeleteRowAndCol() routine
 *      should be compiled.  Note that for this routine to be
 *      compiled, both DELETE and TRANSLATE should be set true.
 *  STRIP
 *      This specifies that the spStripFills() routine should be compiled.
 *  MODIFIED_NODAL
 *      This specifies that the routine that preorders modified node
 *      admittance matrices should be compiled.  This routine results
 *      in greater speed and accuracy if used with this type of
 *      matrix.
 *  QUAD_ELEMENT
 *      This specifies that the routines that allow four related
 *      elements to be entered into the matrix at once should be
 *      compiled.  These elements are usually related to an
 *      admittance.  The routines affected by QUAD_ELEMENT are the
 *      spGetAdmittance, spGetQuad and spGetOnes routines.
 *  TRANSPOSE
 *      This specifies that the routines that solve the matrix as if
 *      it was transposed should be compiled.  These routines are
 *      useful when performing sensitivity analysis using the adjoint
 *      method.
 *  SCALING
 *      This specifies that the routine that performs scaling on the
 *      matrix should be complied.  Scaling is not strongly
 *      supported.  The routine to scale the matrix is provided, but
 *      no routines are provided to scale and descale the RHS and
 *      Solution vectors.  It is suggested that if scaling is desired,
 *      it only be preformed when the pivot order is being chosen [in
 *      spOrderAndFactor()].  This is the only time scaling has
 *      an effect.  The scaling may then either be removed from the
 *      solution by the user or the scaled factors may simply be
 *      thrown away. [NO]
 *  DOCUMENTATION
 *      This specifies that routines that are used to document the
 *      matrix, such as spPrint() and spFileMatrix(), should be
 *      compiled.
 *  DETERMINANT
 *      This specifies that the routine spDeterminant() should be complied.
 *  STABILITY
 *      This specifies that spLargestElement() and spRoundoff() should
 *      be compiled.  These routines are used to check the stability (and
 *      hence the quality of the pivoting) of the factorization by
 *      computing a bound on the size of the element is the matrix E =
 *      A - LU.  If this bound is very high after applying
 *      spOrderAndFactor(), then the pivot threshold should be raised.
 *      If the bound increases greatly after using spFactor(), then the
 *      matrix should probably be reordered.
 *  CONDITION
 *      This specifies that spCondition() and spNorm(), the code that
 *      computes a good estimate of the condition number of the matrix,
 *      should be compiled.
 *  PSEUDOCONDITION
 *      This specifies that spPseudoCondition(), the code that computes
 *      a crude and easily fooled indicator of ill-conditioning in the
 *      matrix, should be compiled.
 *  MULTIPLICATION
 *      This specifies that the routines to multiply the unfactored
 *      matrix by a vector should be compiled.
 *  FORTRAN
 *      This specifies that the FORTRAN interface routines should be
 *      compiled.  When interfacing to FORTRAN programs, the ARRAY_OFFSET
 *      options should be set to NO.
 *  DEBUG
 *      This specifies that additional error checking will be compiled.
 *      The type of error checked are those that are common when the
 *      matrix routines are first integrated into a user's program.  Once
 *      the routines have been integrated in and are running smoothly, this
 *      option should be turned off.
 */

/* Begin options. */
#define  REAL                           YES
#define  EXPANDABLE                     YES
#define  TRANSLATE                      YES
#define  INITIALIZE                     YES
#define  DIAGONAL_PIVOTING              YES
#define  ARRAY_OFFSET                   NOT FORTRAN
#define  MODIFIED_MARKOWITZ             NO
#define  DELETE                         YES
#define  STRIP                          YES
#define  MODIFIED_NODAL                 YES
#define  QUAD_ELEMENT                   YES
#define  TRANSPOSE                      YES
#define  SCALING                        YES
#define  DOCUMENTATION                  YES
#define  MULTIPLICATION                 YES
#define  DETERMINANT                    YES
#define  STABILITY                      YES
#define  CONDITION                      YES
#define  PSEUDOCONDITION                YES
#define  FORTRAN                        NO
#define  DEBUG                          YES

/*
 *  The following options affect Sparse exports and so are exported as a
 *  side effect.  For this reason they use the `sp' prefix.  The boolean
 *  constants YES an NO are not defined in spMatrix.h to avoid conflicts
 *  with user code, so use 0 for NO and 1 for YES.
 */
#endif /* spINSIDE_SPARSE */
#define  spCOMPLEX                      1
#define  spSEPARATED_COMPLEX_VECTORS    0
#ifdef spINSIDE_SPARSE







/*
 *  MATRIX CONSTANTS
 *
 *  These constants are used throughout the sparse matrix routines.  They
 *  should be set to suit the type of matrix being solved.  Recommendations
 *  are given in brackets.
 *
 *  Some terminology should be defined.  The Markowitz row count is the number
 *  of non-zero elements in a row excluding the one being considered as pivot.
 *  There is one Markowitz row count for every row.  The Markowitz column
 *  is defined similarly for columns.  The Markowitz product for an element
 *  is the product of its row and column counts. It is a measure of how much
 *  work would be required on the next step of the factorization if that
 *  element were chosen to be pivot.  A small Markowitz product is desirable.
 *
 *  >>> Constants descriptions:
 *  DEFAULT_THRESHOLD
 *      The relative threshold used if the user enters an invalid
 *      threshold.  Also the threshold used by spFactor() when
 *      calling spOrderAndFactor().  The default threshold should
 *      not be less than or equal to zero nor larger than one. [0.001]
 *  DIAG_PIVOTING_AS_DEFAULT
 *      This indicates whether spOrderAndFactor() should use diagonal
 *      pivoting as default.  This issue only arises when
 *      spOrderAndFactor() is called from spFactor().
 *  SPACE_FOR_ELEMENTS
 *      This number multiplied by the size of the matrix equals the number
 *      of elements for which memory is initially allocated in
 *      spCreate(). [6]
 *  SPACE_FOR_FILL_INS
 *      This number multiplied by the size of the matrix equals the number
 *      of elements for which memory is initially allocated and specifically
 *      reserved for fill-ins in spCreate(). [4]
 *  ELEMENTS_PER_ALLOCATION
 *      The number of matrix elements requested from the malloc utility on
 *      each call to it.  Setting this value greater than 1 reduces the
 *      amount of overhead spent in this system call. On a virtual memory
 *      machine, its good to allocate slightly less than a page worth of
 *      elements at a time (or some multiple thereof).
 *      [For the VAX, for real only use 41, otherwise use 31]
 *  MINIMUM_ALLOCATED_SIZE
 *      The minimum allocated size of a matrix.  Note that this does not
 *      limit the minimum size of a matrix.  This just prevents having to
 *      resize a matrix many times if the matrix is expandable, large and
 *      allocated with an estimated size of zero.  This number should not
 *      be less than one.
 *  EXPANSION_FACTOR
 *      The amount the allocated size of the matrix is increased when it
 *      is expanded.
 *  MAX_MARKOWITZ_TIES
 *      This number is used for two slightly different things, both of which
 *      relate to the search for the best pivot.  First, it is the maximum
 *      number of elements that are Markowitz tied that will be sifted
 *      through when trying to find the one that is numerically the best.
 *      Second, it creates an upper bound on how large a Markowitz product
 *      can be before it eliminates the possibility of early termination
 *      of the pivot search.  In other words, if the product of the smallest
 *      Markowitz product yet found and TIES_MULTIPLIER is greater than
 *      MAX_MARKOWITZ_TIES, then no early termination takes place.
 *      Set MAX_MARKOWITZ_TIES to some small value if no early termination of
 *      the pivot search is desired. An array of RealNumbers is allocated
 *      of size MAX_MARKOWITZ_TIES so it must be positive and shouldn't
 *      be too large.  Active when MODIFIED_MARKOWITZ is 1 (true).  [100]
 *  TIES_MULTIPLIER
 *      Specifies the number of Markowitz ties that are allowed to occur
 *      before the search for the pivot is terminated early.  Set to some
 *      large value if no early termination of the pivot search is desired.
 *      This number is multiplied times the Markowitz product to determine
 *      how many ties are required for early termination.  This means that
 *      more elements will be searched before early termination if a large
 *      number of fill-ins could be created by accepting what is currently
 *      considered the best choice for the pivot.  Active when
 *      MODIFIED_MARKOWITZ is 1 (true).  Setting this number to zero
 *      effectively eliminates all pivoting, which should be avoided.
 *      This number must be positive.  TIES_MULTIPLIER is also used when
 *      diagonal pivoting breaks down. [5]
 *  DEFAULT_PARTITION
 *      Which partition mode is used by spPartition() as default.
 *      Possibilities include
 *          spDIRECT_PARTITION  -- each row used direct addressing, best for
 *              a few relatively dense matrices.
 *          spINDIRECT_PARTITION  -- each row used indirect addressing, best
 *              for a few very sparse matrices.
 *          spAUTO_PARTITION  -- direct or indirect addressing is chosen on
 *              a row-by-row basis, carries a large overhead, but speeds up
 *              both dense and sparse matrices, best if there is a large
 *              number of matrices that can use the same ordering.
 */

/* Begin constants. */
#define  DEFAULT_THRESHOLD              1.0e-3
#define  DIAG_PIVOTING_AS_DEFAULT       YES
#define  SPACE_FOR_ELEMENTS             6
#define  SPACE_FOR_FILL_INS             4
#define  ELEMENTS_PER_ALLOCATION        31
#define  MINIMUM_ALLOCATED_SIZE         6
#define  EXPANSION_FACTOR               1.5
#define  MAX_MARKOWITZ_TIES             100
#define  TIES_MULTIPLIER                5
#define  DEFAULT_PARTITION              spAUTO_PARTITION






/*
 *  PRINTER WIDTH
 *
 *  This macro characterize the printer for the spPrint() routine.
 *
 *  >>> Macros:
 *  PRINTER_WIDTH
 *      The number of characters per page width.  Set to 80 for terminal,
 *      132 for line printer.
 */

/*  Begin printer constants. */
#define  PRINTER_WIDTH  80






/*
 *  MACHINE CONSTANTS
 *
 *  These numbers must be updated when the program is ported to a new machine.
 */

/* Begin machine constants. */

#ifdef notdef  /* __STDC__ */
/*
 * This code is currently deleted because most ANSI standard C compilers
 * do not provide the standard header files yet.
 */
#   include <limits.h>
#   include <float.h>
#   define  MACHINE_RESOLUTION      DBL_EPSILON
#   define  LARGEST_REAL            DBL_MAX
#   define  SMALLEST_REAL           DBL_MIN
#   define  LARGEST_SHORT_INTEGER   SHRT_MAX
#   define  LARGEST_LONG_INTEGER    LONG_MAX
#else /* NOT defined(__STDC__) */

/* VAX machine constants */
#ifdef vax
#   define  MACHINE_RESOLUTION      6.93889e-18
#   define  LARGEST_REAL            1.70141e+38
#   define  SMALLEST_REAL           2.938743e-39
#   define  LARGEST_SHORT_INTEGER   32766
#   define  LARGEST_LONG_INTEGER    2147483646
#endif

/* hp9000 machine constants */
#ifdef hpux
/* These values are correct for hp9000/300.  Should be correct for others. */
#   define  MACHINE_RESOLUTION      8.9e-15
#   define  LARGEST_REAL            1.79769313486231e+308
#   define  SMALLEST_REAL           2.22507385850721e-308
#   define  LARGEST_SHORT_INTEGER   32766
#   define  LARGEST_LONG_INTEGER    2147483646
#endif

/* Sun machine constants */
#ifdef sun
/* These values are rumored to be the correct values. */
#   define  MACHINE_RESOLUTION      8.9e-15
#   define  LARGEST_REAL            1.79769313486231e+308
#   define  SMALLEST_REAL           2.22507385850721e-308
#   define  LARGEST_SHORT_INTEGER   32766
#   define  LARGEST_LONG_INTEGER    2147483646
#endif
#endif /* NOT defined(__STDC__) */






/*
 *  ANNOTATION
 *
 *  This macro changes the amount of annotation produced by the matrix
 *  routines.  The annotation is used as a debugging aid.  Change the number
 *  associated with ANNOTATE to change the amount of annotation produced by
 *  the program.
 */

/* Begin annotation definitions. */
#define  ANNOTATE		NONE

#define  NONE			0
#define  ON_STRANGE_BEHAVIOR	1
#define  FULL			2

#endif /* spINSIDE_SPARSE */
#endif /* spCONFIG_DEFS */
@EOF

chmod 644 spConfig.h

echo x - spDefs.h
cat >spDefs.h <<'@EOF'
/*
 *  DATA STRUCTURE AND MACRO DEFINITIONS for Sparse.
 *
 *  Author:                     Advising professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      UC Berkeley
 *
 *  This file contains common type definitions and macros for the sparse
 *  matrix routines.  These definitions are of no interest to the user.
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 *
 *  $Date: 88/06/18 11:13:40 $
 *  $Revision: 1.2 $
 */




/*
 *  IMPORTS
 */

#include <stdio.h>


/*
 *  If running lint, change some of the compiler options to get a more
 *  complete inspection.
 */

#ifdef lint
#undef  REAL
#undef  spCOMPLEX
#undef  EXPANDABLE
#undef  TRANSLATE
#undef  INITIALIZE
#undef  DELETE
#undef  STRIP
#undef  MODIFIED_NODAL
#undef  QUAD_ELEMENT
#undef  TRANSPOSE
#undef  SCALING
#undef  DOCUMENTATION
#undef  MULTIPLICATION
#undef  DETERMINANT
#undef  CONDITION
#undef  PSEUDOCONDITION
#undef  FORTRAN
#undef  DEBUG

#define  REAL                           YES
#define  spCOMPLEX                      YES
#define  EXPANDABLE                     YES
#define  TRANSLATE                      YES
#define  INITIALIZE                     YES
#define  DELETE                         YES
#define  STRIP                          YES
#define  MODIFIED_NODAL                 YES
#define  QUAD_ELEMENT                   YES
#define  TRANSPOSE                      YES
#define  SCALING                        YES
#define  DOCUMENTATION                  YES
#define  MULTIPLICATION                 YES
#define  DETERMINANT                    YES
#define  CONDITION                      YES
#define  PSEUDOCONDITION                YES
#define  FORTRAN                        YES
#define  DEBUG                          YES

#define  LINT                           YES
#else /* not lint */
#define  LINT                           NO
#endif /* not lint */







/*
 *   MACRO DEFINITIONS
 *
 *   Macros are distinguished by using solely capital letters in their
 *   identifiers.  This contrasts with C defined identifiers which are strictly
 *   lower case, and program variable and procedure names which use both upper
 *   and lower case.
 */

/* Begin macros. */

/* Boolean data type */
#define  BOOLEAN        int
#define  NO             0
#define  YES            1
#define  NOT            !
#define  AND            &&
#define  OR             ||

/* NULL pointer */
#ifndef  NULL
#define  NULL           0
#endif

#define  SPARSE_ID      0x772773        /* Arbitrary (is Sparse on phone). */
#define  IS_SPARSE(matrix)      ((matrix) != NULL &&            \
                                 (matrix)->ID == SPARSE_ID)
#define  IS_VALID(matrix)       ((matrix) != NULL &&            \
                                 (matrix)->ID == SPARSE_ID &&   \
                                 (matrix)->Error >= spOKAY &&   \
                                 (matrix)->Error < spFATAL)
#define  IS_FACTORED(matrix)    ((matrix)->Factored && !(matrix)->NeedsOrdering)

/* Macro commands */
/* Macro functions that return the maximum or minimum independent of type. */
#define  MAX(a,b)           ((a) > (b) ? (a) : (b))
#define  MIN(a,b)           ((a) < (b) ? (a) : (b))

/* Macro function that returns the absolute value of a floating point number. */
#define  ABS(a)             ((a) < 0.0 ? -(a) : (a))

/* Macro function that returns the square of a number. */
#define  SQR(a)             ((a)*(a))

/* Macro procedure that swaps two entities. */
#define  SWAP(type, a, b)   {type swapx; swapx = a; a = b; b = swapx;}

/* Macro function that returns the approx absolute value of a complex number. */
#if spCOMPLEX
#define  ELEMENT_MAG(ptr)   (ABS((ptr)->Real) + ABS((ptr)->Imag))
#else
#define  ELEMENT_MAG(ptr)   ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real)
#endif

/* Complex assignment statements. */
#define  CMPLX_ASSIGN(to,from)  \
{   (to).Real = (from).Real;    \
    (to).Imag = (from).Imag;    \
}
#define  CMPLX_CONJ_ASSIGN(to,from)     \
{   (to).Real = (from).Real;            \
    (to).Imag = -(from).Imag;           \
}
#define  CMPLX_NEGATE_ASSIGN(to,from)   \
{   (to).Real = -(from).Real;           \
    (to).Imag = -(from).Imag;           \
}
#define  CMPLX_CONJ_NEGATE_ASSIGN(to,from)      \
{   (to).Real = -(from).Real;                   \
    (to).Imag = (from).Imag;                    \
}
#define  CMPLX_CONJ(a)  (a).Imag = -(a).Imag
#define  CMPLX_NEGATE(a)        \
{   (a).Real = -(a).Real;       \
    (a).Imag = -(a).Imag;       \
}

/* Macro that returns the approx magnitude (L-1 norm) of a complex number. */
#define  CMPLX_1_NORM(a)        (ABS((a).Real) + ABS((a).Imag))

/* Macro that returns the approx magnitude (L-infinity norm) of a complex. */
#define  CMPLX_INF_NORM(a)      (MAX (ABS((a).Real),ABS((a).Imag)))

/* Macro function that returns the magnitude (L-2 norm) of a complex number. */
#define  CMPLX_2_NORM(a)        (sqrt((a).Real*(a).Real + (a).Imag*(a).Imag))

/* Macro function that performs complex addition. */
#define  CMPLX_ADD(to,from_a,from_b)            \
{   (to).Real = (from_a).Real + (from_b).Real;  \
    (to).Imag = (from_a).Imag + (from_b).Imag;  \
}

/* Macro function that performs complex subtraction. */
#define  CMPLX_SUBT(to,from_a,from_b)           \
{   (to).Real = (from_a).Real - (from_b).Real;  \
    (to).Imag = (from_a).Imag - (from_b).Imag;  \
}

/* Macro function that is equivalent to += operator for complex numbers. */
#define  CMPLX_ADD_ASSIGN(to,from)      \
{   (to).Real += (from).Real;           \
    (to).Imag += (from).Imag;           \
}

/* Macro function that is equivalent to -= operator for complex numbers. */
#define  CMPLX_SUBT_ASSIGN(to,from)     \
{   (to).Real -= (from).Real;           \
    (to).Imag -= (from).Imag;           \
}

/* Macro function that multiplies a complex number by a scalar. */
#define  SCLR_MULT(to,sclr,cmplx)       \
{   (to).Real = (sclr) * (cmplx).Real;  \
    (to).Imag = (sclr) * (cmplx).Imag;  \
}

/* Macro function that multiply-assigns a complex number by a scalar. */
#define  SCLR_MULT_ASSIGN(to,sclr)      \
{   (to).Real *= (sclr);                \
    (to).Imag *= (sclr);                \
}

/* Macro function that multiplies two complex numbers. */
#define  CMPLX_MULT(to,from_a,from_b)           \
{   (to).Real = (from_a).Real * (from_b).Real - \
                (from_a).Imag * (from_b).Imag;  \
    (to).Imag = (from_a).Real * (from_b).Imag + \
                (from_a).Imag * (from_b).Real;  \
}

/* Macro function that implements to *= from for complex numbers. */
#define  CMPLX_MULT_ASSIGN(to,from)             \
{   RealNumber to_real_ = (to).Real;            \
    (to).Real = to_real_ * (from).Real -        \
                (to).Imag * (from).Imag;        \
    (to).Imag = to_real_ * (from).Imag +        \
                (to).Imag * (from).Real;        \
}

/* Macro function that multiplies two complex numbers, the first of which is
 * conjugated. */
#define  CMPLX_CONJ_MULT(to,from_a,from_b)      \
{   (to).Real = (from_a).Real * (from_b).Real + \
                (from_a).Imag * (from_b).Imag;  \
    (to).Imag = (from_a).Real * (from_b).Imag - \
                (from_a).Imag * (from_b).Real;  \
}

/* Macro function that multiplies two complex numbers and then adds them
 * to another. to = add + mult_a * mult_b */
#define  CMPLX_MULT_ADD(to,mult_a,mult_b,add)                   \
{   (to).Real = (mult_a).Real * (mult_b).Real -                 \
                (mult_a).Imag * (mult_b).Imag + (add).Real;     \
    (to).Imag = (mult_a).Real * (mult_b).Imag +                 \
                (mult_a).Imag * (mult_b).Real + (add).Imag;     \
}

/* Macro function that subtracts the product of two complex numbers from
 * another.  to = subt - mult_a * mult_b */
#define  CMPLX_MULT_SUBT(to,mult_a,mult_b,subt)                 \
{   (to).Real = (subt).Real - (mult_a).Real * (mult_b).Real +   \
                              (mult_a).Imag * (mult_b).Imag;    \
    (to).Imag = (subt).Imag - (mult_a).Real * (mult_b).Imag -   \
                              (mult_a).Imag * (mult_b).Real;    \
}

/* Macro function that multiplies two complex numbers and then adds them
 * to another. to = add + mult_a* * mult_b where mult_a* represents mult_a
 * conjugate. */
#define  CMPLX_CONJ_MULT_ADD(to,mult_a,mult_b,add)              \
{   (to).Real = (mult_a).Real * (mult_b).Real +                 \
                (mult_a).Imag * (mult_b).Imag + (add).Real;     \
    (to).Imag = (mult_a).Real * (mult_b).Imag -                 \
                (mult_a).Imag * (mult_b).Real + (add).Imag;     \
}

/* Macro function that multiplies two complex numbers and then adds them
 * to another. to += mult_a * mult_b */
#define  CMPLX_MULT_ADD_ASSIGN(to,from_a,from_b)        \
{   (to).Real += (from_a).Real * (from_b).Real -        \
                 (from_a).Imag * (from_b).Imag;         \
    (to).Imag += (from_a).Real * (from_b).Imag +        \
                 (from_a).Imag * (from_b).Real;         \
}

/* Macro function that multiplies two complex numbers and then subtracts them
 * from another. */
#define  CMPLX_MULT_SUBT_ASSIGN(to,from_a,from_b)       \
{   (to).Real -= (from_a).Real * (from_b).Real -        \
                 (from_a).Imag * (from_b).Imag;         \
    (to).Imag -= (from_a).Real * (from_b).Imag +        \
                 (from_a).Imag * (from_b).Real;         \
}

/* Macro function that multiplies two complex numbers and then adds them
 * to the destination. to += from_a* * from_b where from_a* represents from_a
 * conjugate. */
#define  CMPLX_CONJ_MULT_ADD_ASSIGN(to,from_a,from_b)   \
{   (to).Real += (from_a).Real * (from_b).Real +        \
                 (from_a).Imag * (from_b).Imag;         \
    (to).Imag += (from_a).Real * (from_b).Imag -        \
                 (from_a).Imag * (from_b).Real;         \
}

/* Macro function that multiplies two complex numbers and then subtracts them
 * from the destination. to -= from_a* * from_b where from_a* represents from_a
 * conjugate. */
#define  CMPLX_CONJ_MULT_SUBT_ASSIGN(to,from_a,from_b)  \
{   (to).Real -= (from_a).Real * (from_b).Real +        \
                 (from_a).Imag * (from_b).Imag;         \
    (to).Imag -= (from_a).Real * (from_b).Imag -        \
                 (from_a).Imag * (from_b).Real;         \
}

/*
 * Macro functions that provide complex division.
 */

/* Complex division:  to = num / den */
#define CMPLX_DIV(to,num,den)                                           \
{   RealNumber  r_, s_;                                                 \
    if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR      \
        ((den).Real < (den).Imag AND (den).Real <= -(den).Imag))        \
    {   r_ = (den).Imag / (den).Real;                                   \
        s_ = (den).Real + r_*(den).Imag;                                \
        (to).Real = ((num).Real + r_*(num).Imag)/s_;                    \
        (to).Imag = ((num).Imag - r_*(num).Real)/s_;                    \
    }                                                                   \
    else                                                                \
    {   r_ = (den).Real / (den).Imag;                                   \
        s_ = (den).Imag + r_*(den).Real;                                \
        (to).Real = (r_*(num).Real + (num).Imag)/s_;                    \
        (to).Imag = (r_*(num).Imag - (num).Real)/s_;                    \
    }                                                                   \
}

/* Complex division and assignment:  num /= den */
#define CMPLX_DIV_ASSIGN(num,den)                                       \
{   RealNumber  r_, s_, t_;                                             \
    if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR      \
        ((den).Real < (den).Imag AND (den).Real <= -(den).Imag))        \
    {   r_ = (den).Imag / (den).Real;                                   \
        s_ = (den).Real + r_*(den).Imag;                                \
        t_ = ((num).Real + r_*(num).Imag)/s_;                           \
        (num).Imag = ((num).Imag - r_*(num).Real)/s_;                   \
        (num).Real = t_;                                                \
    }                                                                   \
    else                                                                \
    {   r_ = (den).Real / (den).Imag;                                   \
        s_ = (den).Imag + r_*(den).Real;                                \
        t_ = (r_*(num).Real + (num).Imag)/s_;                           \
        (num).Imag = (r_*(num).Imag - (num).Real)/s_;                   \
        (num).Real = t_;                                                \
    }                                                                   \
}

/* Complex reciprocation:  to = 1.0 / den */
#define CMPLX_RECIPROCAL(to,den)                                        \
{   RealNumber  r_;                                                     \
    if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR      \
        ((den).Real < (den).Imag AND (den).Real <= -(den).Imag))        \
    {   r_ = (den).Imag / (den).Real;                                   \
        (to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \
    }                                                                   \
    else                                                                \
    {   r_ = (den).Real / (den).Imag;                                   \
        (to).Real = -r_*((to).Imag = -1.0/((den).Imag + r_*(den).Real));\
    }                                                                   \
}






/*
 *  ASSERT and ABORT
 *
 *  Macro used to assert that if the code is working correctly, then 
 *  a condition must be true.  If not, then execution is terminated
 *  and an error message is issued stating that there is an internal
 *  error and giving the file and line number.  These assertions are
 *  not evaluated unless the DEBUG flag is true.
 */

#if DEBUG
#define ASSERT(condition) if (NOT(condition)) ABORT()
#else
#define ASSERT(condition)
#endif

#if DEBUG
#define  ABORT()                                                        \
{   (void)fflush(stdout);                                               \
    (void)fprintf(stderr, "sparse: panic in file `%s' at line %d.\n",   \
            __FILE__, __LINE__);                                        \
    (void)fflush(stderr);                                               \
    abort();                                                            \
}
#else
#define  ABORT()
#endif





/*
 *  IMAGINARY VECTORS
 *
 *  The imaginary vectors iRHS and iSolution are only needed when the
 *  options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set.  The following
 *  macro makes it easy to include or exclude these vectors as needed.
 */

#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
#define IMAG_VECTORS    , iRHS, iSolution
#define IMAG_RHS        , iRHS
#else
#define IMAG_VECTORS
#define IMAG_RHS
#endif






/*
 * MEMORY ALLOCATION
 */

extern char *malloc(), *calloc(), *realloc();
#ifdef ultrix
    extern void free();
    extern void abort();
#else
    extern free();
    extern abort();
#endif

#define ALLOC(type,number)  ((type *)malloc((unsigned)(sizeof(type)*(number))))
#define REALLOC(ptr,type,number)  \
           ptr = (type *)realloc((char *)ptr,(unsigned)(sizeof(type)*(number)))
#define FREE(ptr) { if ((ptr) != NULL) free((char *)(ptr)); (ptr) = NULL; }


/* Calloc that properly handles allocating a cleared vector. */
#define CALLOC(ptr,type,number)                         \
{   int i; ptr = ALLOC(type, number);                   \
    if (ptr != (type *)NULL)                            \
        for(i=(number)-1;i>=0; i--) ptr[i] = (type) 0;  \
}







/*
 *  REAL NUMBER
 */

/* Begin `RealNumber'. */

typedef  spREAL  RealNumber, *RealVector;








/*
 *  COMPLEX NUMBER DATA STRUCTURE
 *
 *  >>> Structure fields:
 *  Real  (RealNumber)
 *      The real portion of the number.  Real must be the first
 *      field in this structure.
 *  Imag  (RealNumber)
 *      The imaginary portion of the number. This field must follow
 *      immediately after Real.
 */

/* Begin `ComplexNumber'. */

typedef  struct
{   RealNumber  Real;
    RealNumber  Imag;
} ComplexNumber, *ComplexVector;








/*
 *  MATRIX ELEMENT DATA STRUCTURE
 *
 *  Every nonzero element in the matrix is stored in a dynamically allocated
 *  MatrixElement structure.  These structures are linked together in an
 *  orthogonal linked list.  Two different MatrixElement structures exist.
 *  One is used when only real matrices are expected, it is missing an entry
 *  for imaginary data.  The other is used if complex matrices are expected.
 *  It contains an entry for imaginary data.
 *
 *  >>> Structure fields:
 *  Real  (RealNumber)
 *      The real portion of the value of the element.  Real must be the first
 *      field in this structure.
 *  Imag  (RealNumber)
 *      The imaginary portion of the value of the element. If the matrix
 *      routines are not compiled to handle complex matrices, then this
 *      field does not exist.  If it exists, it must follow immediately after
 *      Real.
 *  Row  (int)
 *      The row number of the element.
 *  Col  (int)
 *      The column number of the element.
 *  NextInRow  (struct MatrixElement *)
 *      NextInRow contains a pointer to the next element in the row to the
 *      right of this element.  If this element is the last nonzero in the
 *      row then NextInRow contains NULL.
 *  NextInCol  (struct MatrixElement *)
 *      NextInCol contains a pointer to the next element in the column below
 *      this element.  If this element is the last nonzero in the column then
 *      NextInCol contains NULL.
 *  pInitInfo  (char *)
 *      Pointer to user data used for initialization of the matrix element.
 *      Initialized to NULL.
 *
 *  >>> Type definitions:
 *  ElementPtr
 *      A pointer to a MatrixElement.
 *  ArrayOfElementPtrs
 *      An array of ElementPtrs.  Used for FirstInRow, FirstInCol and
 *      Diag pointer arrays.
 */

/* Begin `MatrixElement'. */

struct  MatrixElement
{   RealNumber   Real;
#if spCOMPLEX
    RealNumber   Imag;
#endif
    int          Row;
    int          Col;
    struct MatrixElement  *NextInRow;
    struct MatrixElement  *NextInCol;
#if INITIALIZE
    char        *pInitInfo;
#endif
};

typedef  struct MatrixElement  *ElementPtr;
typedef  ElementPtr  *ArrayOfElementPtrs;








/*
 *  ALLOCATION DATA STRUCTURE
 *
 *  The sparse matrix routines keep track of all memory that is allocated by
 *  the operating system so the memory can later be freed.  This is done by
 *  saving the pointers to all the chunks of memory that are allocated to a
 *  particular matrix in an allocation list.  That list is organized as a
 *  linked list so that it can grow without a priori bounds.
 *
 *  >>> Structure fields:
 *  AllocatedPtr  (char *)
 *      Pointer to chunk of memory that has been allocated for the matrix.
 *  NextRecord  (struct  AllocationRecord *)
 *      Pointer to the next allocation record.
 */

/* Begin `AllocationRecord'. */
struct AllocationRecord
{   char  *AllocatedPtr;
    struct  AllocationRecord  *NextRecord;
};

typedef  struct  AllocationRecord  *AllocationListPtr;









/*
 *  FILL-IN LIST DATA STRUCTURE
 *
 *  The sparse matrix routines keep track of all fill-ins separately from
 *  user specified elements so they may be removed by spStripFills().  Fill-ins
 *  are allocated in bunched in what is called a fill-in lists.  The data
 *  structure defined below is used to organize these fill-in lists into a
 *  linked-list.
 *
 *  >>> Structure fields:
 *  pFillinList  (ElementPtr)
 *      Pointer to a fill-in list, or a bunch of fill-ins arranged contiguously
 *      in memory.
 *  NumberOfFillinsInList  (int)
 *      Seems pretty self explanatory to me.
 *  Next  (struct  FillinListNodeStruct *)
 *      Pointer to the next fill-in list structures.
 */

/* Begin `FillinListNodeStruct'. */
struct FillinListNodeStruct
{   ElementPtr  pFillinList;
    int         NumberOfFillinsInList;
    struct      FillinListNodeStruct  *Next;
};










/*
 *  MATRIX FRAME DATA STRUCTURE
 *
 *  This structure contains all the pointers that support the orthogonal
 *  linked list that contains the matrix elements.  Also included in this
 *  structure are other numbers and pointers that are used globally by the
 *  sparse matrix routines and are associated with one particular matrix.
 *
 *  >>> Type definitions:
 *  MatrixPtr
 *      A pointer to MatrixFrame.  Essentially, a pointer to the matrix.
 *
 *  >>> Structure fields:
 *  AbsThreshold  (RealNumber)
 *      The absolute magnitude an element must have to be considered as a
 *      pivot candidate, except as a last resort.
 *  AllocatedExtSize  (int)
 *      The allocated size of the arrays used to translate external row and
 *      column numbers to their internal values.
 *  AllocatedSize  (int)
 *      The currently allocated size of the matrix; the size the matrix can
 *      grow to when EXPANDABLE is set true and AllocatedSize is the largest
 *      the matrix can get without requiring that the matrix frame be
 *      reallocated.
 *  Complex  (BOOLEAN)
 *      The flag which indicates whether the matrix is complex (true) or
 *      real.
 *  CurrentSize  (int)
 *      This number is used during the building of the matrix when the
 *      TRANSLATE option is set true.  It indicates the number of internal
 *      rows and columns that have elements in them.
 *  Diag  (ArrayOfElementPtrs)
 *      Array of pointers that points to the diagonal elements.
 *  DoCmplxDirect  (BOOLEAN *)
 *      Array of flags, one for each column in matrix.  If a flag is true
 *      then corresponding column in a complex matrix should be eliminated
 *      in spFactor() using direct addressing (rather than indirect
 *      addressing).
 *  DoRealDirect  (BOOLEAN *)
 *      Array of flags, one for each column in matrix.  If a flag is true
 *      then corresponding column in a real matrix should be eliminated
 *      in spFactor() using direct addressing (rather than indirect
 *      addressing).
 *  Elements  (int)
 *      The number of original elements (total elements minus fill ins)
 *      present in matrix.
 *  Error  (int)
 *      The error status of the sparse matrix package.
 *  ExtSize  (int)
 *      The value of the largest external row or column number encountered.
 *  ExtToIntColMap  (int [])
 *      An array that is used to convert external columns number to internal
 *      external column numbers.  Present only if TRANSLATE option is set true.
 *  ExtToIntRowMap  (int [])
 *      An array that is used to convert external row numbers to internal
 *      external row numbers.  Present only if TRANSLATE option is set true.
 *  Factored  (BOOLEAN)
 *      Indicates if matrix has been factored.  This flag is set true in
 *      spFactor() and spOrderAndFactor() and set false in spCreate()
 *      and spClear().
 *  Fillins  (int)
 *      The number of fill-ins created during the factorization the matrix.
 *  FirstInCol  (ArrayOfElementPtrs)
 *      Array of pointers that point to the first nonzero element of the
 *      column corresponding to the index.
 *  FirstInRow  (ArrayOfElementPtrs)
 *      Array of pointers that point to the first nonzero element of the row
 *      corresponding to the index.
 *  ID  (unsigned long int)
 *      A constant that provides the sparse data structure with a signature.
 *      When DEBUG is true, all externally available sparse routines check
 *      this signature to assure they are operating on a valid matrix.
 *  Intermediate  (RealVector)
 *      Temporary storage used in the spSolve routines. Intermediate is an
 *      array used during forward and backward substitution.  It is
 *      commonly called y when the forward and backward substitution process is
 *      denoted  Ax = b => Ly = b and Ux = y.
 *  InternalVectorsAllocated  (BOOLEAN)
 *      A flag that indicates whether the markowitz vectors and the
 *      Intermediate vector have been created.
 *      These vectors are created in CreateInternalVectors().
 *  IntToExtColMap  (int [])
 *      An array that is used to convert internal column numbers to external
 *      external column numbers.
 *  IntToExtRowMap  (int [])
 *      An array that is used to convert internal row numbers to external
 *      external row numbers.
 *  MarkowitzCol  (int [])
 *      An array that contains the count of the non-zero elements excluding
 *      the pivots for each column. Used to generate and update MarkowitzProd.
 *  MarkowitzProd  (long [])
 *      The array of the products of the Markowitz row and column counts. The
 *      element with the smallest product is the best pivot to use to maintain
 *      sparsity.
 *  MarkowitzRow  (int [])
 *      An array that contains the count of the non-zero elements excluding
 *      the pivots for each row. Used to generate and update MarkowitzProd.
 *  MaxRowCountInLowerTri  (int)
 *      The maximum number of off-diagonal element in the rows of L, the
 *      lower triangular matrix.  This quantity is used when computing an
 *      estimate of the roundoff error in the matrix.
 *  NeedsOrdering  (BOOLEAN)
 *      This is a flag that signifies that the matrix needs to be ordered
 *      or reordered.  NeedsOrdering is set true in spCreate() and
 *      spGetElement() or spGetAdmittance() if new elements are added to the
 *      matrix after it has been previously factored.  It is set false in
 *      spOrderAndFactor().
 *  NumberOfInterchangesIsOdd  (BOOLEAN)
 *      Flag that indicates the sum of row and column interchange counts
 *      is an odd number.  Used when determining the sign of the determinant.
 *  Partitioned  (BOOLEAN)
 *      This flag indicates that the columns of the matrix have been 
 *      partitioned into two groups.  Those that will be addressed directly
 *      and those that will be addressed indirectly in spFactor().
 *  PivotsOriginalCol  (int)
 *      Column pivot was chosen from.
 *  PivotsOriginalRow  (int)
 *      Row pivot was chosen from.
 *  PivotSelectionMethod  (char)
 *      Character that indicates which pivot search method was successful.
 *  PreviousMatrixWasComplex  (BOOLEAN)
 *      This flag in needed to determine how to clear the matrix.  When
 *      dealing with real matrices, it is important that the imaginary terms
 *      in the matrix elements be zero.  Thus, if the previous matrix was
 *      complex, then the current matrix will be cleared as if it were complex
 *      even if it is real.
 *  RelThreshold  (RealNumber)
 *      The magnitude an element must have relative to others in its row
 *      to be considered as a pivot candidate, except as a last resort.
 *  Reordered  (BOOLEAN)
 *      This flag signifies that the matrix has been reordered.  It
 *      is cleared in spCreate(), set in spMNA_Preorder() and
 *      spOrderAndFactor() and is used in spPrint().
 *  RowsLinked  (BOOLEAN)
 *      A flag that indicates whether the row pointers exist.  The AddByIndex
 *      routines do not generate the row pointers, which are needed by some
 *      of the other routines, such as spOrderAndFactor() and spScale().
 *      The row pointers are generated in the function spcLinkRows().
 *  SingularCol  (int)
 *      Normally zero, but if matrix is found to be singular, SingularCol is
 *      assigned the external column number of pivot that was zero.
 *  SingularRow  (int)
 *      Normally zero, but if matrix is found to be singular, SingularRow is
 *      assigned the external row number of pivot that was zero.
 *  Singletons  (int)
 *      The number of singletons available for pivoting.  Note that if row I
 *      and column I both contain singletons, only one of them is counted.
 *  Size  (int)
 *      Number of rows and columns in the matrix.  Does not change as matrix
 *      is factored.
 *  TrashCan  (MatrixElement)
 *      This is a dummy MatrixElement that is used to by the user to stuff
 *      data related to the zero row or column.  In other words, when the user
 *      adds an element in row zero or column zero, then the matrix returns
 *      a pointer to TrashCan.  In this way the user can have a uniform way
 *      data into the matrix independent of whether a component is connected
 *      to ground.
 *
 *  >>> The remaining fields are related to memory allocation.
 *  TopOfAllocationList  (AllocationListPtr)
 *      Pointer which points to the top entry in a list. The list contains
 *      all the pointers to the segments of memory that have been allocated
 *      to this matrix. This is used when the memory is to be freed on
 *      deallocation of the matrix.
 *  RecordsRemaining  (int)
 *      Number of slots left in the list of allocations.
 *  NextAvailElement  (ElementPtr)
 *      Pointer to the next available element which has been allocated but as
 *      yet is unused. Matrix elements are allocated in groups of
 *      ELEMENTS_PER_ALLOCATION in order to speed element allocation and
 *      freeing.
 *  ElementsRemaining  (int)
 *      Number of unused elements left in last block of elements allocated.
 *  NextAvailFillin  (ElementPtr)
 *      Pointer to the next available fill-in which has been allocated but
 *      as yet is unused.  Fill-ins are allocated in a group in order to keep
 *      them physically close in memory to the rest of the matrix.
 *  FillinsRemaining  (int)
 *      Number of unused fill-ins left in the last block of fill-ins
 *      allocated.
 *  FirstFillinListNode  (FillinListNodeStruct *)
 *      A pointer to the head of the linked-list that keeps track of the
 *      lists of fill-ins.
 *  LastFillinListNode  (FillinListNodeStruct *)
 *      A pointer to the tail of the linked-list that keeps track of the
 *      lists of fill-ins.
 */

/* Begin `MatrixFrame'. */
struct  MatrixFrame
{   RealNumber                   AbsThreshold;
    int                          AllocatedSize;
    int                          AllocatedExtSize;
    BOOLEAN                      Complex;
    int                          CurrentSize;
    ArrayOfElementPtrs           Diag;
    BOOLEAN                     *DoCmplxDirect;
    BOOLEAN                     *DoRealDirect;
    int                          Elements;
    int                          Error;
    int                          ExtSize;
    int                         *ExtToIntColMap;
    int                         *ExtToIntRowMap;
    BOOLEAN                      Factored;
    int                          Fillins;
    ArrayOfElementPtrs           FirstInCol;
    ArrayOfElementPtrs           FirstInRow;
    unsigned long                ID;
    RealVector                   Intermediate;
    BOOLEAN                      InternalVectorsAllocated;
    int                         *IntToExtColMap;
    int                         *IntToExtRowMap;
    int                         *MarkowitzRow;
    int                         *MarkowitzCol;
    long                        *MarkowitzProd;
    int                          MaxRowCountInLowerTri;
    BOOLEAN                      NeedsOrdering;
    BOOLEAN                      NumberOfInterchangesIsOdd;
    BOOLEAN                      Partitioned;
    int                          PivotsOriginalCol;
    int                          PivotsOriginalRow;
    char                         PivotSelectionMethod;
    BOOLEAN                      PreviousMatrixWasComplex;
    RealNumber                   RelThreshold;
    BOOLEAN                      Reordered;
    BOOLEAN                      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;
@EOF

chmod 644 spDefs.h

echo x - spFactor.c
cat >spFactor.c <<'@EOF'
/*
 *  MATRIX FACTORIZATION MODULE
 *
 *  Author:                     Advising Professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      UC Berkeley
 *
 *  This file contains the routines to factor the matrix into LU form.
 *
 *  >>> User accessible functions contained in this file:
 *  spOrderAndFactor
 *  spFactor
 *  spPartition
 *
 *  >>> Other functions contained in this file:
 *  FactorComplexMatrix         CreateInternalVectors
 *  CountMarkowitz              MarkowitzProducts
 *  SearchForPivot              SearchForSingleton
 *  QuicklySearchDiagonal       SearchDiagonal
 *  SearchEntireMatrix          FindLargestInCol
 *  FindBiggestInColExclude     ExchangeRowsAndCols
 *  spcRowExchange              spcColExchange
 *  ExchangeColElements         ExchangeRowElements
 *  RealRowColElimination       ComplexRowColElimination
 *  UpdateMarkowitzNumbers      CreateFillin
 *  MatrixIsSingular            ZeroPivot
 *  WriteStatus
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 */

#ifndef lint
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 $";
#endif



/*
 *  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"






/*
 *  ORDER AND FACTOR MATRIX
 *
 *  This routine chooses a pivot order for the matrix and factors it
 *  into LU form.  It handles both the initial factorization and subsequent
 *  factorizations when a reordering is desired.  This is handled in a manner
 *  that is transparent to the user.  The routine uses a variation of
 *  Gauss's method where the pivots are associated with L and the
 *  diagonal terms of U are one.
 *
 *  >>> Returned:
 *  The error code is returned.  Possible errors are listed below.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *  RHS  <input>  (RealVector)
 *      Representative right-hand side vector that is used to determine
 *      pivoting order when the right hand side vector is sparse.  If
 *      RHS is a NULL pointer then the RHS vector is assumed to
 *      be full and it is not used when determining the pivoting
 *      order.
 *  RelThreshold  <input>  (RealNumber)
 *      This number determines what the pivot relative threshold will
 *      be.  It should be between zero and one.  If it is one then the
 *      pivoting method becomes complete pivoting, which is very slow
 *      and tends to fill up the matrix.  If it is set close to zero
 *      the pivoting method becomes strict Markowitz with no
 *      threshold.  The pivot threshold is used to eliminate pivot
 *      candidates that would cause excessive element growth if they
 *      were used.  Element growth is the cause of roundoff error.
 *      Element growth occurs even in well-conditioned matrices.
 *      Setting the RelThreshold large will reduce element growth and
 *      roundoff error, but setting it too large will cause execution
 *      time to be excessive and will result in a large number of
 *      fill-ins.  If this occurs, accuracy can actually be degraded
 *      because of the large number of operations required on the
 *      matrix due to the large number of fill-ins.  A good value seems
 *      to be 0.001.  The default is chosen by giving a value larger
 *      than one or less than or equal to zero.  This value should be
 *      increased and the matrix resolved if growth is found to be
 *      excessive.  Changing the pivot threshold does not improve
 *      performance on matrices where growth is low, as is often the
 *      case with ill-conditioned matrices.  Once a valid threshold is
 *      given, it becomes the new default.  The default value of
 *      RelThreshold was choosen for use with nearly diagonally
 *      dominant matrices such as node- and modified-node admittance
 *      matrices.  For these matrices it is usually best to use
 *      diagonal pivoting.  For matrices without a strong diagonal, it
 *      is usually best to use a larger threshold, such as 0.01 or
 *      0.1.
 *  AbsThreshold  <input>  (RealNumber)
 *      The absolute magnitude an element must have to be considered
 *      as a pivot candidate, except as a last resort.  This number
 *      should be set significantly smaller than the smallest diagonal
 *      element that is is expected to be placed in the matrix.  If
 *      there is no reasonable prediction for the lower bound on these
 *      elements, then AbsThreshold should be set to zero.
 *      AbsThreshold is used to reduce the possibility of choosing as a
 *      pivot an element that has suffered heavy cancellation and as a
 *      result mainly consists of roundoff error.  Once a valid
 *      threshold is given, it becomes the new default.
 *  DiagPivoting  <input>  (BOOLEAN)
 *      A flag indicating that pivot selection should be confined to the
 *      diagonal if possible.  If DiagPivoting is nonzero and if
 *      DIAGONAL_PIVOTING is enabled pivots will be chosen only from
 *      the diagonal unless there are no diagonal elements that satisfy
 *      the threshold criteria.  Otherwise, the entire reduced
 *      submatrix is searched when looking for a pivot.  The diagonal
 *      pivoting in Sparse is efficient and well refined, while the
 *      off-diagonal pivoting is not.  For symmetric and near symmetric
 *      matrices, it is best to use diagonal pivoting because it
 *      results in the best performance when reordering the matrix and
 *      when factoring the matrix without ordering.  If there is a
 *      considerable amount of nonsymmetry in the matrix, then
 *      off-diagonal pivoting may result in a better equation ordering
 *      simply because there are more pivot candidates to choose from.
 *      A better ordering results in faster subsequent factorizations.
 *      However, the initial pivot selection process takes considerably
 *      longer for off-diagonal pivoting.
 *
 *  >>> Local variables:
 *  pPivot  (ElementPtr)
 *      Pointer to the element being used as a pivot.
 *  ReorderingRequired  (BOOLEAN)
 *      Flag that indicates whether reordering is required.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 *  spSINGULAR
 *  spSMALL_PIVOT
 *  Error is cleared in this function.
 */

int
spOrderAndFactor( eMatrix, RHS, RelThreshold, AbsThreshold, DiagPivoting )

char *eMatrix;
RealNumber  RHS[], RelThreshold, AbsThreshold;
BOOLEAN DiagPivoting;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
ElementPtr  pPivot;
int  Step, Size, ReorderingRequired;
ElementPtr SearchForPivot();
RealNumber LargestInCol, FindLargestInCol();

/* Begin `spOrderAndFactor'. */
    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);

    Matrix->Error = spOKAY;
    Size = Matrix->Size;
    if (RelThreshold <= 0.0) RelThreshold = Matrix->RelThreshold;
    if (RelThreshold > 1.0) RelThreshold = Matrix->RelThreshold;
    Matrix->RelThreshold = RelThreshold;
    if (AbsThreshold < 0.0) AbsThreshold = Matrix->AbsThreshold;
    Matrix->AbsThreshold = AbsThreshold;
    ReorderingRequired = NO;

    if (NOT Matrix->NeedsOrdering)
    {
/* Matrix has been factored before and reordering is not required. */
        for (Step = 1; Step <= Size; Step++)
        {   pPivot = Matrix->Diag[Step];
            LargestInCol = FindLargestInCol(pPivot->NextInCol);
            if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot)))
            {   if (Matrix->Complex)
                    ComplexRowColElimination( Matrix, pPivot );
                else
                    RealRowColElimination( Matrix, pPivot );
            }
            else
            {   ReorderingRequired = YES;
                break; /* for loop */
            }
        }
        if (NOT ReorderingRequired)
            goto Done;
        else
        {
/*
 * A pivot was not large enough to maintain accuracy,
 * so a partial reordering is required.
 */

#if ANNOTATE >= ON_STRANGE_BEHAVIOR
            printf("Reordering,  Step = %1d\n", Step);
#endif
        }
    } /* End of if(NOT Matrix->NeedsOrdering) */
    else
    {
/*
 * This is the first time the matrix has been factored.  These few statements
 * indicate to the rest of the code that a full reodering is required rather
 * than a partial reordering, which occurs during a failure of a fast
 * factorization.
 */
        Step = 1;
        if (NOT Matrix->RowsLinked)
            spcLinkRows( Matrix );
        if (NOT Matrix->InternalVectorsAllocated)
            CreateInternalVectors( Matrix );
        if (Matrix->Error >= spFATAL)
            return Matrix->Error;
    }

/* Form initial Markowitz products. */
    CountMarkowitz( Matrix, RHS, Step );
    MarkowitzProducts( Matrix, Step );
    Matrix->MaxRowCountInLowerTri = -1;

/* Perform reordering and factorization. */
    for (; Step <= Size; Step++)
    {   pPivot = SearchForPivot( Matrix, Step, DiagPivoting );
        if (pPivot == NULL) return MatrixIsSingular( Matrix, Step );
        ExchangeRowsAndCols( Matrix, pPivot, Step );

        if (Matrix->Complex)
            ComplexRowColElimination( Matrix, pPivot );
        else
            RealRowColElimination( Matrix, pPivot );

        if (Matrix->Error >= spFATAL) return Matrix->Error;
        UpdateMarkowitzNumbers( Matrix, pPivot );

#if ANNOTATE == FULL
        WriteStatus( Matrix, Step );
#endif
    }

Done:
    Matrix->NeedsOrdering = NO;
    Matrix->Reordered = YES;
    Matrix->Factored = YES;

    return Matrix->Error;
}







/*
 *  FACTOR MATRIX
 *
 *  This routine is the companion routine to spOrderAndFactor().
 *  Unlike spOrderAndFactor(), spFactor() cannot change the ordering.
 *  It is also faster than spOrderAndFactor().  The standard way of
 *  using these two routines is to first use spOrderAndFactor() for the
 *  initial factorization.  For subsequent factorizations, spFactor()
 *  is used if there is some assurance that little growth will occur
 *  (say for example, that the matrix is diagonally dominant).  If
 *  spFactor() is called for the initial factorization of the matrix,
 *  then spOrderAndFactor() is automatically called with the default
 *  threshold.  This routine uses "row at a time" LU factorization.
 *  Pivots are associated with the lower triangular matrix and the
 *  diagonals of the upper triangular matrix are ones.
 *
 *  >>> Returned:
 *  The error code is returned.  Possible errors are listed below.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 *  spSINGULAR
 *  spZERO_DIAG
 *  spSMALL_PIVOT
 *  Error is cleared in this function.
 */

int
spFactor( eMatrix )

char *eMatrix;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement;
register  ElementPtr  pColumn;
register  int  Step, Size;
RealNumber Mult;

/* Begin `spFactor'. */
    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);

    if (Matrix->NeedsOrdering)
    {   return spOrderAndFactor( eMatrix, (RealVector)NULL,
                                 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT );
    }
    if (NOT Matrix->Partitioned) spPartition( eMatrix, spDEFAULT_PARTITION );
#if spCOMPLEX
    if (Matrix->Complex) return FactorComplexMatrix( Matrix );
#endif

#if REAL
    Size = Matrix->Size;

    if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 );
    Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real;

/* Start factorization. */
    for (Step = 2; Step <= Size; Step++)
    {   if (Matrix->DoRealDirect[Step])
        {   /* Update column using direct addressing scatter-gather. */
            register RealNumber *Dest = (RealNumber *)Matrix->Intermediate;

/* Scatter. */
            pElement = Matrix->FirstInCol[Step];
            while (pElement != NULL)
            {   Dest[pElement->Row] = pElement->Real;
                pElement = pElement->NextInCol;
            }

/* Update column. */
            pColumn = Matrix->FirstInCol[Step];
            while (pColumn->Row < Step)
            {   pElement = Matrix->Diag[pColumn->Row];
                pColumn->Real = Dest[pColumn->Row] * pElement->Real;
                while ((pElement = pElement->NextInCol) != NULL)
                    Dest[pElement->Row] -= pColumn->Real * pElement->Real;
                pColumn = pColumn->NextInCol;
            }

/* Gather. */
            pElement = Matrix->Diag[Step]->NextInCol;
            while (pElement != NULL)
            {   pElement->Real = Dest[pElement->Row];
                pElement = pElement->NextInCol;
            }

/* Check for singular matrix. */
            if (Dest[Step] == 0.0) return ZeroPivot( Matrix, Step );
            Matrix->Diag[Step]->Real = 1.0 / Dest[Step];
        }
        else
        {   /* Update column using indirect addressing scatter-gather. */
            register RealNumber **pDest = (RealNumber **)Matrix->Intermediate;

/* Scatter. */
            pElement = Matrix->FirstInCol[Step];
            while (pElement != NULL)
            {   pDest[pElement->Row] = &pElement->Real;
                pElement = pElement->NextInCol;
            }

/* Update column. */
            pColumn = Matrix->FirstInCol[Step];
            while (pColumn->Row < Step)
            {   pElement = Matrix->Diag[pColumn->Row];
                Mult = (*pDest[pColumn->Row] *= pElement->Real);
                while ((pElement = pElement->NextInCol) != NULL)
                    *pDest[pElement->Row] -= Mult * pElement->Real;
                pColumn = pColumn->NextInCol;
            }

/* Check for singular matrix. */
            if (Matrix->Diag[Step]->Real == 0.0)
                return ZeroPivot( Matrix, Step );
            Matrix->Diag[Step]->Real = 1.0 / Matrix->Diag[Step]->Real;
        }
    }

    Matrix->Factored = YES;
    return (Matrix->Error = spOKAY);
#endif /* REAL */
}






#if spCOMPLEX
/*
 *  FACTOR COMPLEX MATRIX
 *
 *  This routine is the companion routine to spFactor(), it
 *  handles complex matrices.  It is otherwise identical.
 *
 *  >>> Returned:
 *  The error code is returned.  Possible errors are listed below.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *
 *  >>> Possible errors:
 *  spSINGULAR
 *  Error is cleared in this function.
 */

static int
FactorComplexMatrix( Matrix )

MatrixPtr  Matrix;
{
register  ElementPtr  pElement;
register  ElementPtr  pColumn;
register  int  Step, Size;
ComplexNumber Mult, Pivot;

/* Begin `FactorComplexMatrix'. */
    ASSERT(Matrix->Complex);

    Size = Matrix->Size;
    pElement = Matrix->Diag[1];
    if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, 1 );
/* Cmplx expr: *pPivot = 1.0 / *pPivot. */
    CMPLX_RECIPROCAL( *pElement, *pElement );

/* Start factorization. */
    for (Step = 2; Step <= Size; Step++)
    {   if (Matrix->DoCmplxDirect[Step])
        {   /* Update column using direct addressing scatter-gather. */
            register  ComplexNumber  *Dest;
            Dest = (ComplexNumber *)Matrix->Intermediate;

/* Scatter. */
            pElement = Matrix->FirstInCol[Step];
            while (pElement != NULL)
            {   Dest[pElement->Row] = *(ComplexNumber *)pElement;
                pElement = pElement->NextInCol;
            }

/* Update column. */
            pColumn = Matrix->FirstInCol[Step];
            while (pColumn->Row < Step)
            {   pElement = Matrix->Diag[pColumn->Row];
                /* Cmplx expr: Mult = Dest[pColumn->Row] * (1.0 / *pPivot). */
                CMPLX_MULT(Mult, Dest[pColumn->Row], *pElement);
                CMPLX_ASSIGN(*pColumn, Mult);
                while ((pElement = pElement->NextInCol) != NULL)
                {   /* Cmplx expr: Dest[pElement->Row] -= Mult * pElement */
                    CMPLX_MULT_SUBT_ASSIGN(Dest[pElement->Row],Mult,*pElement);
                }
                pColumn = pColumn->NextInCol;
            }

/* Gather. */
            pElement = Matrix->Diag[Step]->NextInCol;
            while (pElement != NULL)
            {   *(ComplexNumber *)pElement = Dest[pElement->Row];
                pElement = pElement->NextInCol;
            }

/* Check for singular matrix. */
            Pivot = Dest[Step];
            if (CMPLX_1_NORM(Pivot) == 0.0) return ZeroPivot( Matrix, Step );
            CMPLX_RECIPROCAL( *Matrix->Diag[Step], Pivot );  
        }
        else
        {   /* Update column using direct addressing scatter-gather. */
            register  ComplexNumber  **pDest;
            pDest = (ComplexNumber **)Matrix->Intermediate;

/* Scatter. */
            pElement = Matrix->FirstInCol[Step];
            while (pElement != NULL)
            {   pDest[pElement->Row] = (ComplexNumber *)pElement;
                pElement = pElement->NextInCol;
            }

/* Update column. */
            pColumn = Matrix->FirstInCol[Step];
            while (pColumn->Row < Step)
            {   pElement = Matrix->Diag[pColumn->Row];
                /* Cmplx expr: Mult = *pDest[pColumn->Row] * (1.0 / *pPivot). */
                CMPLX_MULT(Mult, *pDest[pColumn->Row], *pElement);
                CMPLX_ASSIGN(*pDest[pColumn->Row], Mult);
                while ((pElement = pElement->NextInCol) != NULL)
                {  /* Cmplx expr: *pDest[pElement->Row] -= Mult * pElement */
                   CMPLX_MULT_SUBT_ASSIGN(*pDest[pElement->Row],Mult,*pElement);
                }
                pColumn = pColumn->NextInCol;
            }

/* Check for singular matrix. */
            pElement = Matrix->Diag[Step];
            if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, Step );
            CMPLX_RECIPROCAL( *pElement, *pElement );  
        }
    }

    Matrix->Factored = YES;
    return (Matrix->Error = spOKAY);
}
#endif /* spCOMPLEX */






/*
 *  PARTITION MATRIX
 *
 *  This routine determines the cost to factor each row using both
 *  direct and indirect addressing and decides, on a row-by-row basis,
 *  which addressing mode is fastest.  This information is used in
 *  spFactor() to speed the factorization.
 *
 *  When factoring a previously ordered matrix using spFactor(), Sparse
 *  operates on a row-at-a-time basis.  For speed, on each step, the
 *  row being updated is copied into a full vector and the operations
 *  are performed on that vector.  This can be done one of two ways,
 *  either using direct addressing or indirect addressing.  Direct
 *  addressing is fastest when the matrix is relatively dense and
 *  indirect addressing is best when the matrix is quite sparse.  The
 *  user selects the type of partition used with Mode.  If Mode is set
 *  to spDIRECT_PARTITION, then the all rows are placed in the direct
 *  addressing partition.  Similarly, if Mode is set to
 *  spINDIRECT_PARTITION, then the all rows are placed in the indirect
 *  addressing partition.  By setting Mode to spAUTO_PARTITION, the
 *  user allows Sparse to select the partition for each row
 *  individually.  spFactor() generally runs faster if Sparse is
 *  allowed to choose its own partitioning, however choosing a
 *  partition is expensive.  The time required to choose a partition is
 *  of the same order of the cost to factor the matrix.  If you plan to
 *  factor a large number of matrices with the same structure, it is
 *  best to let Sparse choose the partition.  Otherwise, you should
 *  choose the partition based on the predicted density of the matrix.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *  Mode  <input>  (int)
 *      Mode must be one of three special codes: spDIRECT_PARTITION,
 *      spINDIRECT_PARTITION, or spAUTO_PARTITION.
 */

void
spPartition( eMatrix, Mode )

char *eMatrix;
int Mode;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement, pColumn;
register  int  Step, Size;
register  int  *Nc, *No, *Nm;
BOOLEAN *DoRealDirect, *DoCmplxDirect;

/* Begin `spPartition'. */
    ASSERT( IS_SPARSE( Matrix ) );
    if (Matrix->Partitioned) return;
    Size = Matrix->Size;
    DoRealDirect = Matrix->DoRealDirect;
    DoCmplxDirect = Matrix->DoCmplxDirect;
    Matrix->Partitioned = YES;

/* If partition is specified by the user, this is easy. */
    if (Mode == spDEFAULT_PARTITION) Mode = DEFAULT_PARTITION;
    if (Mode == spDIRECT_PARTITION)
    {   for (Step = 1; Step <= Size; Step++)
#if REAL
            DoRealDirect[Step] = YES;
#endif
#if spCOMPLEX
            DoCmplxDirect[Step] = YES;
#endif
        return;
    }
    else if (Mode == spINDIRECT_PARTITION)
    {   for (Step = 1; Step <= Size; Step++)
#if REAL
            DoRealDirect[Step] = NO;
#endif
#if spCOMPLEX
            DoCmplxDirect[Step] = NO;
#endif
        return;
    }
    else ASSERT( Mode == spAUTO_PARTITION );

/* Otherwise, count all operations needed in when factoring matrix. */
    Nc = (int *)Matrix->MarkowitzRow;
    No = (int *)Matrix->MarkowitzCol;
    Nm = (int *)Matrix->MarkowitzProd;

/* Start mock-factorization. */
    for (Step = 1; Step <= Size; Step++)
    {   Nc[Step] = No[Step] = Nm[Step] = 0;

        pElement = Matrix->FirstInCol[Step];
        while (pElement != NULL)
        {   Nc[Step]++;
            pElement = pElement->NextInCol;
        }

        pColumn = Matrix->FirstInCol[Step];
        while (pColumn->Row < Step)
        {   pElement = Matrix->Diag[pColumn->Row];
            Nm[Step]++;
            while ((pElement = pElement->NextInCol) != NULL)
                No[Step]++;
            pColumn = pColumn->NextInCol;
        }
    }

    for (Step = 1; Step <= Size; Step++)
    {
/*
 * The following are just estimates based on a count on the number of
 * machine instructions used on each machine to perform the various
 * tasks.  It was assumed that each machine instruction required the
 * same amount of time (I don't believe this is true for the VAX, and
 * have no idea if this is true for the 68000 family).  For optimum
 * performance, these numbers should be tuned to the machine.
 *   Nc is the number of nonzero elements in the column.
 *   Nm is the number of multipliers in the column.
 *   No is the number of operations in the inner loop.
 */

#define generic
#ifdef hp9000s300
#if REAL
        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
#endif
#if spCOMPLEX
        /* On the hp350, it is never profitable to use direct for complex. */
        DoCmplxDirect[Step] = NO;
#endif
#undef generic
#endif

#ifdef vax
#if REAL
        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
#endif
#if spCOMPLEX
        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
#endif
#undef generic
#endif

#ifdef generic
#if REAL
        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
#endif
#if spCOMPLEX
        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
#endif
#undef generic
#endif
    }

#if (ANNOTATE == FULL)
    {   int Ops = 0;
        for (Step = 1; Step <= Size; Step++)
            Ops += No[Step];
        printf("Operation count for inner loop of factorization = %d.\n", Ops);
    }
#endif
    return;
}







/*
 *  CREATE INTERNAL VECTORS
 *
 *  Creates the Markowitz and Intermediate vectors.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *
 *  >>> Local variables:
 *  SizePlusOne  (unsigned)
 *      Size of the arrays to be allocated.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static
CreateInternalVectors( Matrix )

MatrixPtr  Matrix;
{
int  Size;

/* Begin `CreateInternalVectors'. */
/* Create Markowitz arrays. */
    Size= Matrix->Size;

    if (Matrix->MarkowitzRow == NULL)
    {   if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL)
            Matrix->Error = spNO_MEMORY;
    }
    if (Matrix->MarkowitzCol == NULL)
    {   if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL)
            Matrix->Error = spNO_MEMORY;
    }
    if (Matrix->MarkowitzProd == NULL)
    {   if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL)
            Matrix->Error = spNO_MEMORY;
    }

/* Create DoDirect vectors for use in spFactor(). */
#if REAL
    if (Matrix->DoRealDirect == NULL)
    {   if (( Matrix->DoRealDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
            Matrix->Error = spNO_MEMORY;
    }
#endif
#if spCOMPLEX
    if (Matrix->DoCmplxDirect == NULL)
    {   if (( Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
            Matrix->Error = spNO_MEMORY;
    }
#endif

/* Create Intermediate vectors for use in MatrixSolve. */
#if spCOMPLEX
    if (Matrix->Intermediate == NULL)
    {   if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL)
            Matrix->Error = spNO_MEMORY;
    }
#else
    if (Matrix->Intermediate == NULL)
    {   if ((Matrix->Intermediate = ALLOC(RealNumber, Size+1)) == NULL)
            Matrix->Error = spNO_MEMORY;
    }
#endif

    if (Matrix->Error != spNO_MEMORY)
        Matrix->InternalVectorsAllocated = YES;
    return;
}







/*
 *  COUNT MARKOWITZ
 *
 *  Scans Matrix to determine the Markowitz counts for each row and column.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  RHS  <input>  (RealVector)
 *      Representative right-hand side vector that is used to determine
 *      pivoting order when the right hand side vector is sparse.  If
 *      RHS is a NULL pointer then the RHS vector is assumed to be full
 *      and it is not used when determining the pivoting order.
 *  Step  <input>  (int)
 *     Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  Count  (int)
 *     Temporary counting variable.
 *  ExtRow  (int)
 *     The external row number that corresponds to I.
 *  pElement  (ElementPtr)
 *     Pointer to matrix elements.
 *  Size  (int)
 *     The size of the matrix.
 */

static
CountMarkowitz( Matrix, RHS, Step )

MatrixPtr Matrix;
register RealVector  RHS;
int Step;
{
register int  Count, I, Size = Matrix->Size;
register ElementPtr  pElement;
int  ExtRow;

/* Begin `CountMarkowitz'. */

/* Correct array pointer for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX
        if (RHS != NULL) --RHS;
#else
        if (RHS != NULL)
        {   if (Matrix->Complex) RHS -= 2;
            else --RHS;
        }
#endif
#endif

/* Generate MarkowitzRow Count for each row. */
    for (I = Step; I <= Size; I++)
    {
/* Set Count to -1 initially to remove count due to pivot element. */
        Count = -1;
        pElement = Matrix->FirstInRow[I];
        while (pElement != NULL AND pElement->Col < Step)
            pElement = pElement->NextInRow;
        while (pElement != NULL)
        {   Count++;
            pElement = pElement->NextInRow;
        }

/* Include nonzero elements in the RHS vector. */
        ExtRow = Matrix->IntToExtRowMap[I];

#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX
        if (RHS != NULL)
            if (RHS[ExtRow] != 0.0)  Count++;
#else
        if (RHS != NULL)
        {   if (Matrix->Complex)
            {   if ((RHS[2*ExtRow] != 0.0) OR (RHS[2*ExtRow+1] != 0.0))
                    Count++;
            }
            else if (RHS[I] != 0.0) Count++;
        }
#endif
        Matrix->MarkowitzRow[I] = Count;
    }

/* Generate the MarkowitzCol count for each column. */
    for (I = Step; I <= Size; I++)
    {
/* Set Count to -1 initially to remove count due to pivot element. */
        Count = -1;
        pElement = Matrix->FirstInCol[I];
        while (pElement != NULL AND pElement->Row < Step)
            pElement = pElement->NextInCol;
        while (pElement != NULL)
        {   Count++;
            pElement = pElement->NextInCol;
        }
        Matrix->MarkowitzCol[I] = Count;
    }
    return;
}










/*
 *  MARKOWITZ PRODUCTS
 *
 *  Calculates MarkowitzProduct for each diagonal element from the Markowitz
 *  counts.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local Variables:
 *  pMarkowitzProduct  (long *)
 *      Pointer that points into MarkowitzProduct array. Is used to
 *      sequentially access entries quickly.
 *  pMarkowitzRow  (int *)
 *      Pointer that points into MarkowitzRow array. Is used to sequentially
 *      access entries quickly.
 *  pMarkowitzCol  (int *)
 *      Pointer that points into MarkowitzCol array. Is used to sequentially
 *      access entries quickly.
 *  Product  (long)
 *      Temporary storage for Markowitz product./
 *  Size  (int)
 *      The size of the matrix.
 */

static
MarkowitzProducts( Matrix, Step )

MatrixPtr Matrix;
int Step;
{
register  int  I, *pMarkowitzRow, *pMarkowitzCol;
register  long  Product, *pMarkowitzProduct;
register  int  Size = Matrix->Size;
double fProduct;

/* Begin `MarkowitzProducts'. */
    Matrix->Singletons = 0;

    pMarkowitzProduct = &(Matrix->MarkowitzProd[Step]);
    pMarkowitzRow = &(Matrix->MarkowitzRow[Step]);
    pMarkowitzCol = &(Matrix->MarkowitzCol[Step]);

    for (I = Step; I <= Size; I++)
    {
/* If chance of overflow, use real numbers. */
        if ((*pMarkowitzRow > LARGEST_SHORT_INTEGER AND *pMarkowitzCol != 0) OR
            (*pMarkowitzCol > LARGEST_SHORT_INTEGER AND *pMarkowitzRow != 0))
        {   fProduct = (double)(*pMarkowitzRow++) * (double)(*pMarkowitzCol++);
            if (fProduct >= LARGEST_LONG_INTEGER)
                *pMarkowitzProduct++ = LARGEST_LONG_INTEGER;
            else
                *pMarkowitzProduct++ = fProduct;
        }
        else
        {   Product = *pMarkowitzRow++ * *pMarkowitzCol++;
            if ((*pMarkowitzProduct++ = Product) == 0)
                Matrix->Singletons++;
        }
    }
    return;
}











/*
 *  SEARCH FOR BEST PIVOT
 *
 *  Performs a search to determine the element with the lowest Markowitz
 *  Product that is also acceptable.  An acceptable element is one that is
 *  larger than the AbsThreshold and at least as large as RelThreshold times
 *  the largest element in the same column.  The first step is to look for
 *  singletons if any exist.  If none are found, then all the diagonals are
 *  searched. The diagonal is searched once quickly using the assumption that
 *  elements on the diagonal are large compared to other elements in their
 *  column, and so the pivot can be chosen only on the basis of the Markowitz
 *  criterion.  After a element has been chosen to be pivot on the basis of
 *  its Markowitz product, it is checked to see if it is large enough.
 *  Waiting to the end of the Markowitz search to check the size of a pivot
 *  candidate saves considerable time, but is not guaranteed to find an
 *  acceptable pivot.  Thus if unsuccessful a second pass of the diagonal is
 *  made.  This second pass checks to see if an element is large enough during
 *  the search, not after it.  If still no acceptable pivot candidate has
 *  been found, the search expands to cover the entire matrix.
 *
 *  >>> Returned:
 *  A pointer to the element chosen to be pivot.  If every element in the
 *  matrix is zero, then NULL is returned.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      The row and column number of the beginning of the reduced submatrix.
 *
 *  >>> Local variables:
 *  ChosenPivot  (ElementPtr)
 *      Pointer to element that has been chosen to be the pivot.
 *
 *  >>> Possible errors:
 *  spSINGULAR
 *  spSMALL_PIVOT
 */

static ElementPtr
SearchForPivot( Matrix, Step, DiagPivoting )

MatrixPtr Matrix;
int Step, DiagPivoting;
{
register ElementPtr  ChosenPivot;
ElementPtr  SearchForSingleton();
ElementPtr  QuicklySearchDiagonal();
ElementPtr  SearchDiagonal();
ElementPtr  SearchEntireMatrix();

/* Begin `SearchForPivot'. */

/* If singletons exist, look for an acceptable one to use as pivot. */
    if (Matrix->Singletons)
    {   ChosenPivot = SearchForSingleton( Matrix, Step );
        if (ChosenPivot != NULL)
        {   Matrix->PivotSelectionMethod = 's';
            return ChosenPivot;
        }
    }

#if DIAGONAL_PIVOTING
    if (DiagPivoting)
    {
/*
 * Either no singletons exist or they weren't acceptable.  Take quick first
 * pass at searching diagonal.  First search for element on diagonal of 
 * remaining submatrix with smallest Markowitz product, then check to see
 * if it okay numerically.  If not, QuicklySearchDiagonal fails.
 */
        ChosenPivot = QuicklySearchDiagonal( Matrix, Step );
        if (ChosenPivot != NULL)
        {   Matrix->PivotSelectionMethod = 'q';
            return ChosenPivot;
        }

/*
 * Quick search of diagonal failed, carefully search diagonal and check each
 * pivot candidate numerically before even tentatively accepting it.
 */
        ChosenPivot = SearchDiagonal( Matrix, Step );
        if (ChosenPivot != NULL)
        {   Matrix->PivotSelectionMethod = 'd';
            return ChosenPivot;
        }
    }
#endif /* DIAGONAL_PIVOTING */

/* No acceptable pivot found yet, search entire matrix. */
    ChosenPivot = SearchEntireMatrix( Matrix, Step );
    Matrix->PivotSelectionMethod = 'e';

    return ChosenPivot;
}









/*
 *  SEARCH FOR SINGLETON TO USE AS PIVOT
 *
 *  Performs a search to find a singleton to use as the pivot.  The
 *  first acceptable singleton is used.  A singleton is acceptable if
 *  it is larger in magnitude than the AbsThreshold and larger
 *  than RelThreshold times the largest of any other elements in the same
 *  column.  It may seem that a singleton need not satisfy the
 *  relative threshold criterion, however it is necessary to prevent
 *  excessive growth in the RHS from resulting in overflow during the
 *  forward and backward substitution.  A singleton does not need to
 *  be on the diagonal to be selected.
 *
 *  >>> Returned:
 *  A pointer to the singleton chosen to be pivot.  In no singleton is
 *  acceptable, return NULL.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  ChosenPivot  (ElementPtr)
 *      Pointer to element that has been chosen to be the pivot.
 *  PivotMag  (RealNumber)
 *      Magnitude of ChosenPivot.
 *  Singletons  (int)
 *      The count of the number of singletons that can be used as pivots.
 *      A local version of Matrix->Singletons.
 *  pMarkowitzProduct  (long *)
 *      Pointer that points into MarkowitzProduct array. It is used to quickly
 *      access successive Markowitz products.
 */

static ElementPtr
SearchForSingleton( Matrix, Step )

MatrixPtr Matrix;
int Step;
{
register  ElementPtr  ChosenPivot;
register  int  I;
register  long  *pMarkowitzProduct;
int  Singletons;
RealNumber  PivotMag, FindBiggestInColExclude();

/* Begin `SearchForSingleton'. */
/* Initialize pointer that is to scan through MarkowitzProduct vector. */
    pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+1]);
    Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];

/* Decrement the count of available singletons, on the assumption that an
 * acceptable one will be found. */
    Singletons = Matrix->Singletons--;

/*
 * Assure that following while loop will always terminate, this is just
 * preventive medicine, if things are working right this should never
 * be needed.
 */
    Matrix->MarkowitzProd[Step-1] = 0;

    while (Singletons-- > 0)
    {
/* Singletons exist, find them. */

/*
 * This is tricky.  Am using a pointer to sequentially step through the
 * MarkowitzProduct array.  Search terminates when singleton (Product = 0)
 * is found.  Note that the conditional in the while statement
 * ( *pMarkowitzProduct ) is true as long as the MarkowitzProduct is not
 * equal to zero.  The row (and column) index on the diagonal is then
 * calculated by subtracting the pointer to the Markowitz product of
 * the first diagonal from the pointer to the Markowitz product of the
 * desired element, the singleton.
 *
 * Search proceeds from the end (high row and column numbers) to the
 * beginning (low row and column numbers) so that rows and columns with
 * large Markowitz products will tend to be move to the bottom of the
 * matrix.  However, choosing Diag[Step] is desirable because it would
 * require no row and column interchanges, so inspect it first by
 * putting its Markowitz product at the end of the MarkowitzProd
 * vector.
 */

        while ( *pMarkowitzProduct-- )
        {   /*
             * N bottles of beer on the wall;
             * N bottles of beer.
             * you take one down and pass it around;
             * N-1 bottles of beer on the wall.
             */
        }
        I = pMarkowitzProduct - Matrix->MarkowitzProd + 1;

/* Assure that I is valid. */
        if (I < Step) break;  /* while (Singletons-- > 0) */
        if (I > Matrix->Size) I = Step;

/* Singleton has been found in either/both row or/and column I. */
        if ((ChosenPivot = Matrix->Diag[I]) != NULL)
        {
/* Singleton lies on the diagonal. */
            PivotMag = ELEMENT_MAG(ChosenPivot);
            if
            (    PivotMag > Matrix->AbsThreshold AND
                 PivotMag > Matrix->RelThreshold *
                            FindBiggestInColExclude( Matrix, ChosenPivot, Step )
            ) return ChosenPivot;
        }
        else
        {
/* Singleton does not lie on diagonal, find it. */
            if (Matrix->MarkowitzCol[I] == 0)
            {   ChosenPivot = Matrix->FirstInCol[I];
                while ((ChosenPivot != NULL) AND (ChosenPivot->Row < Step))
                    ChosenPivot = ChosenPivot->NextInCol;
                PivotMag = ELEMENT_MAG(ChosenPivot);
                if
                (    PivotMag > Matrix->AbsThreshold AND
                     PivotMag > Matrix->RelThreshold *
                                FindBiggestInColExclude( Matrix, ChosenPivot,
                                                         Step )
                ) return ChosenPivot;
                else
                {   if (Matrix->MarkowitzRow[I] == 0)
                    {   ChosenPivot = Matrix->FirstInRow[I];
                        while((ChosenPivot != NULL) AND (ChosenPivot->Col<Step))
                            ChosenPivot = ChosenPivot->NextInRow;
                        PivotMag = ELEMENT_MAG(ChosenPivot);
                        if
                        (    PivotMag > Matrix->AbsThreshold AND
                             PivotMag > Matrix->RelThreshold *
                                        FindBiggestInColExclude( Matrix,
                                                                 ChosenPivot,
                                                                 Step )
                        ) return ChosenPivot;
                    }
                }
            }
            else
            {   ChosenPivot = Matrix->FirstInRow[I];
                while ((ChosenPivot != NULL) AND (ChosenPivot->Col < Step))
                    ChosenPivot = ChosenPivot->NextInRow;
                PivotMag = ELEMENT_MAG(ChosenPivot);
                if
                (    PivotMag > Matrix->AbsThreshold AND
                     PivotMag > Matrix->RelThreshold *
                                FindBiggestInColExclude( Matrix, ChosenPivot,
                                                         Step )
                ) return ChosenPivot;
            }
        }
/* Singleton not acceptable (too small), try another. */
    } /* end of while(lSingletons>0) */

/*
 * All singletons were unacceptable.  Restore Matrix->Singletons count.
 * Initial assumption that an acceptable singleton would be found was wrong.
 */
    Matrix->Singletons++;
    return NULL;
}












#if DIAGONAL_PIVOTING
#if MODIFIED_MARKOWITZ
/*
 *  QUICK SEARCH OF DIAGONAL FOR PIVOT WITH MODIFIED MARKOWITZ CRITERION
 *
 *  Searches the diagonal looking for the best pivot.  For a pivot to be
 *  acceptable it must be larger than the pivot RelThreshold times the largest
 *  element in its reduced column.  Among the acceptable diagonals, the
 *  one with the smallest MarkowitzProduct is sought.  Search terminates
 *  early if a diagonal is found with a MarkowitzProduct of one and its
 *  magnitude is larger than the other elements in its row and column.
 *  Since its MarkowitzProduct is one, there is only one other element in
 *  both its row and column, and, as a condition for early termination,
 *  these elements must be located symmetricly in the matrix.  If a tie
 *  occurs between elements of equal MarkowitzProduct, then the element with
 *  the largest ratio between its magnitude and the largest element in its
 *  column is used.  The search will be terminated after a given number of
 *  ties have occurred and the best (largest ratio) of the tied element will
 *  be used as the pivot.  The number of ties that will trigger an early
 *  termination is MinMarkowitzProduct * TIES_MULTIPLIER.
 *
 *  >>> Returned:
 *  A pointer to the diagonal element chosen to be pivot.  If no diagonal is
 *  acceptable, a NULL is returned.
 *
 *  >>> Arguments:
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  ChosenPivot  (ElementPtr)
 *      Pointer to the element that has been chosen to be the pivot.
 *  LargestOffDiagonal  (RealNumber)
 *      Magnitude of the largest of the off-diagonal terms associated with
 *      a diagonal with MarkowitzProduct equal to one.
 *  Magnitude  (RealNumber)
 *      Absolute value of diagonal element.
 *  MaxRatio  (RealNumber)
 *      Among the elements tied with the smallest Markowitz product, MaxRatio
 *      is the best (smallest) ratio of LargestInCol to the diagonal Magnitude
 *      found so far.  The smaller the ratio, the better numerically the
 *      element will be as pivot.
 *  MinMarkowitzProduct  (long)
 *      Smallest Markowitz product found of pivot candidates that lie along
 *      diagonal.
 *  NumberOfTies  (int)
 *      A count of the number of Markowitz ties that have occurred at current
 *      MarkowitzProduct.
 *  pDiag  (ElementPtr)
 *      Pointer to current diagonal element.
 *  pMarkowitzProduct  (long *)
 *      Pointer that points into MarkowitzProduct array. It is used to quickly
 *      access successive Markowitz products.
 *  Ratio  (RealNumber)
 *      For the current pivot candidate, Ratio is the ratio of the largest
 *      element in its column (excluding itself) to its magnitude.
 *  TiedElements  (ElementPtr[])
 *      Array of pointers to the elements with the minimum Markowitz
 *      product.
 *  pOtherInCol  (ElementPtr)
 *      When there is only one other element in a column other than the
 *      diagonal, pOtherInCol is used to point to it.  Used when Markowitz
 *      product is to determine if off diagonals are placed symmetricly.
 *  pOtherInRow  (ElementPtr)
 *      When there is only one other element in a row other than the diagonal,
 *      pOtherInRow is used to point to it.  Used when Markowitz product is
 *      to determine if off diagonals are placed symmetricly.
 */

static ElementPtr
QuicklySearchDiagonal( Matrix, Step )

MatrixPtr Matrix;
int Step;
{
register long  MinMarkowitzProduct, *pMarkowitzProduct;
register  ElementPtr  pDiag, pOtherInRow, pOtherInCol;
int  I, NumberOfTies;
ElementPtr  ChosenPivot, TiedElements[MAX_MARKOWITZ_TIES + 1];
RealNumber  Magnitude, LargestInCol, Ratio, MaxRatio;
RealNumber  LargestOffDiagonal;
RealNumber  FindBiggestInColExclude();

/* Begin `QuicklySearchDiagonal'. */
    NumberOfTies = -1;
    MinMarkowitzProduct = LARGEST_LONG_INTEGER;
    pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]);
    Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];

/* Assure that following while loop will always terminate. */
    Matrix->MarkowitzProd[Step-1] = -1;

/*
 * This is tricky.  Am using a pointer in the inner while loop to
 * sequentially step through the MarkowitzProduct array.  Search
 * terminates when the Markowitz product of zero placed at location
 * Step-1 is found.  The row (and column) index on the diagonal is then
 * calculated by subtracting the pointer to the Markowitz product of
 * the first diagonal from the pointer to the Markowitz product of the
 * desired element. The outer for loop is infinite, broken by using
 * break.
 *
 * Search proceeds from the end (high row and column numbers) to the
 * beginning (low row and column numbers) so that rows and columns with
 * large Markowitz products will tend to be move to the bottom of the
 * matrix.  However, choosing Diag[Step] is desirable because it would
 * require no row and column interchanges, so inspect it first by
 * putting its Markowitz product at the end of the MarkowitzProd
 * vector.
 */

    for(;;)  /* Endless for loop. */
    {   while (MinMarkowitzProduct < *(--pMarkowitzProduct))
        {   /*
             * N bottles of beer on the wall;
             * N bottles of beer.
             * You take one down and pass it around;
             * N-1 bottles of beer on the wall.
             */
        }

        I = pMarkowitzProduct - Matrix->MarkowitzProd;

/* Assure that I is valid; if I < Step, terminate search. */
        if (I < Step) break; /* Endless for loop */
        if (I > Matrix->Size) I = Step;

        if ((pDiag = Matrix->Diag[I]) == NULL)
            continue; /* Endless for loop */
        if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
            continue; /* Endless for loop */

        if (*pMarkowitzProduct == 1)
        {
/* Case where only one element exists in row and column other than diagonal. */

/* Find off diagonal elements. */
            pOtherInRow = pDiag->NextInRow;
            pOtherInCol = pDiag->NextInCol;
            if (pOtherInRow == NULL AND pOtherInCol == NULL)
            {    pOtherInRow = Matrix->FirstInRow[I];
                 while(pOtherInRow != NULL)
                 {   if (pOtherInRow->Col >= Step AND pOtherInRow->Col != I)
                         break;
                     pOtherInRow = pOtherInRow->NextInRow;
                 }
                 pOtherInCol = Matrix->FirstInCol[I];
                 while(pOtherInCol != NULL)
                 {   if (pOtherInCol->Row >= Step AND pOtherInCol->Row != I)
                         break;
                     pOtherInCol = pOtherInCol->NextInCol;
                 }
            }

/* Accept diagonal as pivot if diagonal is larger than off diagonals and the
 * off diagonals are placed symmetricly. */
            if (pOtherInRow != NULL  AND  pOtherInCol != NULL)
            {   if (pOtherInRow->Col == pOtherInCol->Row)
                {   LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow),
                                                      ELEMENT_MAG(pOtherInCol));
                    if (Magnitude >= LargestOffDiagonal)
                    {
/* Accept pivot, it is unlikely to contribute excess error. */
                        return pDiag;
                    }
                }
            }
        }

        if (*pMarkowitzProduct < MinMarkowitzProduct)
        {
/* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */
            TiedElements[0] = pDiag;
            MinMarkowitzProduct = *pMarkowitzProduct;
            NumberOfTies = 0;
        }
        else
        {
/* This case handles Markowitz ties. */
            if (NumberOfTies < MAX_MARKOWITZ_TIES)
            {   TiedElements[++NumberOfTies] = pDiag;
                if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
                    break; /* Endless for loop */
            }
        }
    } /* End of endless for loop. */

/* Test to see if any element was chosen as a pivot candidate. */
    if (NumberOfTies < 0)
        return NULL;

/* Determine which of tied elements is best numerically. */
    ChosenPivot = NULL;
    MaxRatio = 1.0 / Matrix->RelThreshold;

    for (I = 0; I <= NumberOfTies; I++)
    {   pDiag = TiedElements[I];
        Magnitude = ELEMENT_MAG(pDiag);
        LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step );
        Ratio = LargestInCol / Magnitude;
        if (Ratio < MaxRatio)
        {   ChosenPivot = pDiag;
            MaxRatio = Ratio;
        }
    }
    return ChosenPivot;
}










#else /* Not MODIFIED_MARKOWITZ */
/*
 *  QUICK SEARCH OF DIAGONAL FOR PIVOT WITH CONVENTIONAL MARKOWITZ
 *  CRITERION
 *
 *  Searches the diagonal looking for the best pivot.  For a pivot to be
 *  acceptable it must be larger than the pivot RelThreshold times the largest
 *  element in its reduced column.  Among the acceptable diagonals, the
 *  one with the smallest MarkowitzProduct is sought.  Search terminates
 *  early if a diagonal is found with a MarkowitzProduct of one and its
 *  magnitude is larger than the other elements in its row and column.
 *  Since its MarkowitzProduct is one, there is only one other element in
 *  both its row and column, and, as a condition for early termination,
 *  these elements must be located symmetricly in the matrix.
 *
 *  >>> Returned:
 *  A pointer to the diagonal element chosen to be pivot.  If no diagonal is
 *  acceptable, a NULL is returned.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  ChosenPivot  (ElementPtr)
 *      Pointer to the element that has been chosen to be the pivot.
 *  LargestOffDiagonal  (RealNumber)
 *      Magnitude of the largest of the off-diagonal terms associated with
 *      a diagonal with MarkowitzProduct equal to one.
 *  Magnitude  (RealNumber)
 *      Absolute value of diagonal element.
 *  MinMarkowitzProduct  (long)
 *      Smallest Markowitz product found of pivot candidates which are
 *      acceptable.
 *  pDiag  (ElementPtr)
 *      Pointer to current diagonal element.
 *  pMarkowitzProduct  (long *)
 *      Pointer that points into MarkowitzProduct array. It is used to quickly
 *      access successive Markowitz products.
 *  pOtherInCol  (ElementPtr)
 *      When there is only one other element in a column other than the
 *      diagonal, pOtherInCol is used to point to it.  Used when Markowitz
 *      product is to determine if off diagonals are placed symmetricly.
 *  pOtherInRow  (ElementPtr)
 *      When there is only one other element in a row other than the diagonal,
 *      pOtherInRow is used to point to it.  Used when Markowitz product is
 *      to determine if off diagonals are placed symmetricly.
 */

static ElementPtr
QuicklySearchDiagonal( Matrix, Step )

MatrixPtr Matrix;
int Step;
{
register long  MinMarkowitzProduct, *pMarkowitzProduct;
register  ElementPtr  pDiag;
int  I;
ElementPtr  ChosenPivot, pOtherInRow, pOtherInCol;
RealNumber  Magnitude, LargestInCol, LargestOffDiagonal;
RealNumber  FindBiggestInColExclude();

/* Begin `QuicklySearchDiagonal'. */
    ChosenPivot = NULL;
    MinMarkowitzProduct = LARGEST_LONG_INTEGER;
    pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]);
    Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step];

/* Assure that following while loop will always terminate. */
    Matrix->MarkowitzProd[Step-1] = -1;

/*
 * This is tricky.  Am using a pointer in the inner while loop to
 * sequentially step through the MarkowitzProduct array.  Search
 * terminates when the Markowitz product of zero placed at location
 * Step-1 is found.  The row (and column) index on the diagonal is then
 * calculated by subtracting the pointer to the Markowitz product of
 * the first diagonal from the pointer to the Markowitz product of the
 * desired element. The outer for loop is infinite, broken by using
 * break.
 *
 * Search proceeds from the end (high row and column numbers) to the
 * beginning (low row and column numbers) so that rows and columns with
 * large Markowitz products will tend to be move to the bottom of the
 * matrix.  However, choosing Diag[Step] is desirable because it would
 * require no row and column interchanges, so inspect it first by
 * putting its Markowitz product at the end of the MarkowitzProd
 * vector.
 */

    for (;;)  /* Endless for loop. */
    {   while (*(--pMarkowitzProduct) >= MinMarkowitzProduct)
        {   /* Just passing through. */
        }

        I = pMarkowitzProduct - Matrix->MarkowitzProd;

/* Assure that I is valid; if I < Step, terminate search. */
        if (I < Step) break; /* Endless for loop */
        if (I > Matrix->Size) I = Step;

        if ((pDiag = Matrix->Diag[I]) == NULL)
            continue; /* Endless for loop */
        if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
            continue; /* Endless for loop */

        if (*pMarkowitzProduct == 1)
        {
/* Case where only one element exists in row and column other than diagonal. */

/* Find off-diagonal elements. */
            pOtherInRow = pDiag->NextInRow;
            pOtherInCol = pDiag->NextInCol;
            if (pOtherInRow == NULL AND pOtherInCol == NULL)
            {    pOtherInRow = Matrix->FirstInRow[I];
                 while(pOtherInRow != NULL)
                 {   if (pOtherInRow->Col >= Step AND pOtherInRow->Col != I)
                         break;
                     pOtherInRow = pOtherInRow->NextInRow;
                 }
                 pOtherInCol = Matrix->FirstInCol[I];
                 while(pOtherInCol != NULL)
                 {   if (pOtherInCol->Row >= Step AND pOtherInCol->Row != I)
                         break;
                     pOtherInCol = pOtherInCol->NextInCol;
                 }
            }

/* Accept diagonal as pivot if diagonal is larger than off-diagonals and the
 * off-diagonals are placed symmetricly. */
            if (pOtherInRow != NULL  AND  pOtherInCol != NULL)
            {   if (pOtherInRow->Col == pOtherInCol->Row)
                {   LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow),
                                                      ELEMENT_MAG(pOtherInCol));
                    if (Magnitude >= LargestOffDiagonal)
                    {
/* Accept pivot, it is unlikely to contribute excess error. */
                        return pDiag;
                    }
                }
            }
        }

        MinMarkowitzProduct = *pMarkowitzProduct;
        ChosenPivot = pDiag;
    }  /* End of endless for loop. */

    if (ChosenPivot != NULL)
    {   LargestInCol = FindBiggestInColExclude( Matrix, ChosenPivot, Step );
        if( ELEMENT_MAG(ChosenPivot) <= Matrix->RelThreshold * LargestInCol )
            ChosenPivot = NULL;
    }
    return ChosenPivot;
}
#endif /* Not MODIFIED_MARKOWITZ */









/*
 *  SEARCH DIAGONAL FOR PIVOT WITH MODIFIED MARKOWITZ CRITERION
 *
 *  Searches the diagonal looking for the best pivot.  For a pivot to be
 *  acceptable it must be larger than the pivot RelThreshold times the largest
 *  element in its reduced column.  Among the acceptable diagonals, the
 *  one with the smallest MarkowitzProduct is sought.  If a tie occurs
 *  between elements of equal MarkowitzProduct, then the element with
 *  the largest ratio between its magnitude and the largest element in its
 *  column is used.  The search will be terminated after a given number of
 *  ties have occurred and the best (smallest ratio) of the tied element will
 *  be used as the pivot.  The number of ties that will trigger an early
 *  termination is MinMarkowitzProduct * TIES_MULTIPLIER.
 *
 *  >>> Returned:
 *  A pointer to the diagonal element chosen to be pivot.  If no diagonal is
 *  acceptable, a NULL is returned.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  ChosenPivot  (ElementPtr)
 *      Pointer to the element that has been chosen to be the pivot.
 *  Size  (int)
 *      Local version of size which is placed in a register to increase speed.
 *  Magnitude  (RealNumber)
 *      Absolute value of diagonal element.
 *  MinMarkowitzProduct  (long)
 *      Smallest Markowitz product found of those pivot candidates which are
 *      acceptable.
 *  NumberOfTies  (int)
 *      A count of the number of Markowitz ties that have occurred at current
 *      MarkowitzProduct.
 *  pDiag  (ElementPtr)
 *      Pointer to current diagonal element.
 *  pMarkowitzProduct  (long*)
 *      Pointer that points into MarkowitzProduct array. It is used to quickly
 *      access successive Markowitz products.
 *  Ratio  (RealNumber)
 *      For the current pivot candidate, Ratio is the
 *      Ratio of the largest element in its column to its magnitude.
 *  RatioOfAccepted  (RealNumber)
 *      For the best pivot candidate found so far, RatioOfAccepted is the
 *      Ratio of the largest element in its column to its magnitude.
 */

static ElementPtr
SearchDiagonal( Matrix, Step )

MatrixPtr Matrix;
register int Step;
{
register  int  J;
register long  MinMarkowitzProduct, *pMarkowitzProduct;
register  int  I;
register  ElementPtr  pDiag;
int  NumberOfTies, Size = Matrix->Size;
ElementPtr  ChosenPivot;
RealNumber  Magnitude, Ratio, RatioOfAccepted, LargestInCol;
RealNumber  FindBiggestInColExclude();

/* Begin `SearchDiagonal'. */
    ChosenPivot = NULL;
    MinMarkowitzProduct = LARGEST_LONG_INTEGER;
    pMarkowitzProduct = &(Matrix->MarkowitzProd[Size+2]);
    Matrix->MarkowitzProd[Size+1] = Matrix->MarkowitzProd[Step];

/* Start search of diagonal. */
    for (J = Size+1; J > Step; J--)
    {
        if (*(--pMarkowitzProduct) > MinMarkowitzProduct)
            continue; /* for loop */
        if (J > Matrix->Size)
            I = Step;
        else
            I = J;
        if ((pDiag = Matrix->Diag[I]) == NULL)
            continue; /* for loop */
        if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold)
            continue; /* for loop */

/* Test to see if diagonal's magnitude is acceptable. */
        LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step );
        if (Magnitude <= Matrix->RelThreshold * LargestInCol)
            continue; /* for loop */

        if (*pMarkowitzProduct < MinMarkowitzProduct)
        {
/* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */
            ChosenPivot = pDiag;
            MinMarkowitzProduct = *pMarkowitzProduct;
            RatioOfAccepted = LargestInCol / Magnitude;
            NumberOfTies = 0;
        }
        else
        {
/* This case handles Markowitz ties. */
            NumberOfTies++;
            Ratio = LargestInCol / Magnitude;
            if (Ratio < RatioOfAccepted)
            {   ChosenPivot = pDiag;
                RatioOfAccepted = Ratio;
            }
            if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
                return ChosenPivot;
        }
    } /* End of for(Step) */
    return ChosenPivot;
}
#endif /* DIAGONAL_PIVOTING */










/*
 *  SEARCH ENTIRE MATRIX FOR BEST PIVOT
 *
 *  Performs a search over the entire matrix looking for the acceptable
 *  element with the lowest MarkowitzProduct.  If there are several that
 *  are tied for the smallest MarkowitzProduct, the tie is broken by using
 *  the ratio of the magnitude of the element being considered to the largest
 *  element in the same column.  If no element is acceptable then the largest
 *  element in the reduced submatrix is used as the pivot and the
 *  matrix is declared to be spSMALL_PIVOT.  If the largest element is
 *  zero, the matrix is declared to be spSINGULAR.
 *
 *  >>> Returned:
 *  A pointer to the diagonal element chosen to be pivot.  If no element is
 *  found, then NULL is returned and the matrix is spSINGULAR.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  ChosenPivot  (ElementPtr)
 *      Pointer to the element that has been chosen to be the pivot.
 *  LargestElementMag  (RealNumber)
 *      Magnitude of the largest element yet found in the reduced submatrix.
 *  Size  (int)
 *      Local version of Size; placed in a register for speed.
 *  Magnitude  (RealNumber)
 *      Absolute value of diagonal element.
 *  MinMarkowitzProduct  (long)
 *      Smallest Markowitz product found of pivot candidates which are
 *      acceptable.
 *  NumberOfTies  (int)
 *      A count of the number of Markowitz ties that have occurred at current
 *      MarkowitzProduct.
 *  pElement  (ElementPtr)
 *      Pointer to current element.
 *  pLargestElement  (ElementPtr)
 *      Pointer to the largest element yet found in the reduced submatrix.
 *  Product  (long)
 *      Markowitz product for the current row and column.
 *  Ratio  (RealNumber)
 *      For the current pivot candidate, Ratio is the
 *      Ratio of the largest element in its column to its magnitude.
 *  RatioOfAccepted  (RealNumber)
 *      For the best pivot candidate found so far, RatioOfAccepted is the
 *      Ratio of the largest element in its column to its magnitude.
 *
 *  >>> Possible errors:
 *  spSINGULAR
 *  spSMALL_PIVOT
 */

static ElementPtr
SearchEntireMatrix( Matrix, Step )

MatrixPtr Matrix;
int Step;
{
register  int  I, Size = Matrix->Size;
register  ElementPtr  pElement;
int  NumberOfTies;
long  Product, MinMarkowitzProduct;
ElementPtr  ChosenPivot, pLargestElement;
RealNumber  Magnitude, LargestElementMag, Ratio, RatioOfAccepted, LargestInCol;
RealNumber  FindLargestInCol();

/* Begin `SearchEntireMatrix'. */
    ChosenPivot = NULL;
    LargestElementMag = 0.0;
    MinMarkowitzProduct = LARGEST_LONG_INTEGER;

/* Start search of matrix on column by column basis. */
    for (I = Step; I <= Size; I++)
    {   pElement = Matrix->FirstInCol[I];

        while (pElement != NULL AND pElement->Row < Step)
            pElement = pElement->NextInCol;

        if((LargestInCol = FindLargestInCol(pElement)) == 0.0)
            continue; /* for loop */

        while (pElement != NULL)
        {
/* Check to see if element is the largest encountered so far.  If so, record
   its magnitude and address. */
            if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElementMag)
            {   LargestElementMag = Magnitude;
                pLargestElement = pElement;
            }
/* Calculate element's MarkowitzProduct. */
            Product = Matrix->MarkowitzRow[pElement->Row] *
                      Matrix->MarkowitzCol[pElement->Col];

/* Test to see if element is acceptable as a pivot candidate. */
            if ((Product <= MinMarkowitzProduct) AND
                (Magnitude > Matrix->RelThreshold * LargestInCol) AND
                (Magnitude > Matrix->AbsThreshold))
            {
/* Test to see if element has lowest MarkowitzProduct yet found, or whether it
   is tied with an element found earlier. */
                if (Product < MinMarkowitzProduct)
                {
/* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */
                    ChosenPivot = pElement;
                    MinMarkowitzProduct = Product;
                    RatioOfAccepted = LargestInCol / Magnitude;
                    NumberOfTies = 0;
                }
                else
                {
/* This case handles Markowitz ties. */
                    NumberOfTies++;
                    Ratio = LargestInCol / Magnitude;
                    if (Ratio < RatioOfAccepted)
                    {   ChosenPivot = pElement;
                        RatioOfAccepted = Ratio;
                    }
                    if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER)
                        return ChosenPivot;
                }
            }
            pElement = pElement->NextInCol;
        }  /* End of while(pElement != NULL) */
    } /* End of for(Step) */

    if (ChosenPivot != NULL) return ChosenPivot;

    if (LargestElementMag == 0.0)
    {   Matrix->Error = spSINGULAR;
        return NULL;
    }

    Matrix->Error = spSMALL_PIVOT;
    return pLargestElement;
}












/*
 *  DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
 *
 *  This routine searches a column and returns the magnitude of the largest
 *  element.  This routine begins the search at the element pointed to by
 *  pElement, the parameter.
 *
 *  The search is conducted by starting at the element specified by a pointer,
 *  which should be one below the diagonal, and moving down the column.  On
 *  the way down the column, the magnitudes of the elements are tested to see
 *  if they are the largest yet found.
 *
 *  >>> Returned:
 *  The magnitude of the largest element in the column below and including
 *  the one pointed to by the input parameter.
 *
 *  >>> Arguments:
 *  pElement  <input>  (ElementPtr)
 *      The pointer to the first element to be tested.  Also, used by the
 *      routine to access all lower elements in the column.
 *
 *  >>> Local variables:
 *  Largest  (RealNumber)
 *      The magnitude of the largest element.
 *  Magnitude  (RealNumber)
 *      The magnitude of the currently active element.
 */

static RealNumber
FindLargestInCol( pElement )

register  ElementPtr  pElement;
{
RealNumber  Magnitude, Largest = 0.0;

/* Begin `FindLargestInCol'. */
/* Search column for largest element beginning at Element. */
    while (pElement != NULL)
    {   if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
            Largest = Magnitude;
        pElement = pElement->NextInCol;
    }

    return Largest;
}










/*
 *  DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
 *  EXCLUDING AN ELEMENT
 *
 *  This routine searches a column and returns the magnitude of the largest
 *  element.  One given element is specifically excluded from the search.
 *
 *  The search is conducted by starting at the first element in the column
 *  and moving down the column until the active part of the matrix is entered,
 *  i.e. the reduced submatrix.  The rest of the column is then traversed
 *  looking for the largest element.
 *
 *  >>> Returned:
 *  The magnitude of the largest element in the active portion of the column,
 *  excluding the specified element, is returned.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  pElement  <input>  (ElementPtr)
 *      The pointer to the element that is to be excluded from search. Column
 *      to be searched is one that contains this element.  Also used to
 *      access the elements in the column.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.  Indicates where
 *      the active part of the matrix begins.
 *
 *  >>> Local variables:
 *  Col  (int)
 *      The number of the column to be searched.  Also the column number of
 *      the element to be avoided in the search.
 *  Largest  (RealNumber)
 *      The magnitude of the largest element.
 *  Magnitude  (RealNumber)
 *      The magnitude of the currently active element.
 *  Row  (int)
 *      The row number of element to be excluded from the search.
 */

static RealNumber
FindBiggestInColExclude( Matrix, pElement, Step )

MatrixPtr Matrix;
register  ElementPtr  pElement;
register  int Step;
{
register  int  Row;
int  Col;
RealNumber  Largest, Magnitude;

/* Begin `FindBiggestInColExclude'. */
    Row = pElement->Row;
    Col = pElement->Col;
    pElement = Matrix->FirstInCol[Col];

/* Travel down column until reduced submatrix is entered. */
    while ((pElement != NULL) AND (pElement->Row < Step))
        pElement = pElement->NextInCol;

/* Initialize the variable Largest. */
    if (pElement->Row != Row)
        Largest = ELEMENT_MAG(pElement);
    else
        Largest = 0.0;

/* Search rest of column for largest element, avoiding excluded element. */
    while ((pElement = pElement->NextInCol) != NULL)
    {   if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
        {   if (pElement->Row != Row)
                Largest = Magnitude;
        }
    }

    return Largest;
}










/*
 *  EXCHANGE ROWS AND COLUMNS
 *
 *  Exchanges two rows and two columns so that the selected pivot is moved to
 *  the upper left corner of the remaining submatrix.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  pPivot  <input>  (ElementPtr)
 *      Pointer to the current pivot.
 *  Step  <input>  (int)
 *      Index of the diagonal currently being eliminated.
 *
 *  >>> Local variables:
 *  Col  (int)
 *      Column where the pivot was found.
 *  Row  (int)
 *      Row where the pivot was found.
 *  OldMarkowitzProd_Col  (long)
 *      Markowitz product associated with the diagonal element in the row
 *      the pivot was found in.
 *  OldMarkowitzProd_Row  (long)
 *      Markowitz product associated with the diagonal element in the column
 *      the pivot was found in.
 *  OldMarkowitzProd_Step  (long)
 *      Markowitz product associated with the diagonal element that is being
 *      moved so that the pivot can be placed in the upper left-hand corner
 *      of the reduced submatrix.
 */

static
ExchangeRowsAndCols( Matrix, pPivot, Step )

MatrixPtr Matrix;
ElementPtr  pPivot;
register int Step;
{
register  int   Row, Col;
long  OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col;
ElementPtr spcFindElementInCol();

/* Begin `ExchangeRowsAndCols'. */
    Row = pPivot->Row;
    Col = pPivot->Col;
    Matrix->PivotsOriginalRow = Row;
    Matrix->PivotsOriginalCol = Col;

    if ((Row == Step) AND (Col == Step)) return;

/* Exchange rows and columns. */
    if (Row == Col)
    {   spcRowExchange( Matrix, Step, Row );
        spcColExchange( Matrix, Step, Col );
        SWAP( long, Matrix->MarkowitzProd[Step], Matrix->MarkowitzProd[Row] );
        SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Step] );
    }
    else
    {

/* Initialize variables that hold old Markowitz products. */
        OldMarkowitzProd_Step = Matrix->MarkowitzProd[Step];
        OldMarkowitzProd_Row = Matrix->MarkowitzProd[Row];
        OldMarkowitzProd_Col = Matrix->MarkowitzProd[Col];

/* Exchange rows. */
        if (Row != Step)
        {   spcRowExchange( Matrix, Step, Row );
            Matrix->NumberOfInterchangesIsOdd =
                                       NOT Matrix->NumberOfInterchangesIsOdd;
            Matrix->MarkowitzProd[Row] = Matrix->MarkowitzRow[Row] *
                                                   Matrix->MarkowitzCol[Row];

/* Update singleton count. */
            if ((Matrix->MarkowitzProd[Row]==0) != (OldMarkowitzProd_Row==0))
            {   if (OldMarkowitzProd_Row == 0)
                    Matrix->Singletons--;
                else
                    Matrix->Singletons++;
            }
        }

/* Exchange columns. */
        if (Col != Step)
        {   spcColExchange( Matrix, Step, Col );
            Matrix->NumberOfInterchangesIsOdd =
                                       NOT Matrix->NumberOfInterchangesIsOdd;
            Matrix->MarkowitzProd[Col] = Matrix->MarkowitzCol[Col] *
                                                   Matrix->MarkowitzRow[Col];

/* Update singleton count. */
            if ((Matrix->MarkowitzProd[Col]==0) != (OldMarkowitzProd_Col==0))
            {   if (OldMarkowitzProd_Col == 0)
                    Matrix->Singletons--;
                else
                    Matrix->Singletons++;
            }

            Matrix->Diag[Col] = spcFindElementInCol( Matrix,
                                                     Matrix->FirstInCol+Col,
                                                     Col, Col, NO );
        }
        if (Row != Step)
        {   Matrix->Diag[Row] = spcFindElementInCol( Matrix,
                                                     Matrix->FirstInCol+Row,
                                                     Row, Row, NO );
        }
        Matrix->Diag[Step] = spcFindElementInCol( Matrix,
                                                  Matrix->FirstInCol+Step,
                                                  Step, Step, NO );

/* Update singleton count. */
        Matrix->MarkowitzProd[Step] = Matrix->MarkowitzCol[Step] *
                                                    Matrix->MarkowitzRow[Step];
        if ((Matrix->MarkowitzProd[Step]==0) != (OldMarkowitzProd_Step==0))
        {   if (OldMarkowitzProd_Step == 0)
                Matrix->Singletons--;
            else
                Matrix->Singletons++;
        }
    }
    return;
}









/*
 *  EXCHANGE ROWS
 *
 *  Performs all required operations to exchange two rows. Those operations
 *  include: swap FirstInRow pointers, fixing up the NextInCol pointers,
 *  swapping row indexes in MatrixElements, and swapping Markowitz row
 *  counts.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  Row1  <input>  (int)
 *      Row index of one of the rows, becomes the smallest index.
 *  Row2  <input>  (int)
 *      Row index of the other row, becomes the largest index.
 *
 *  Local variables:
 *  Column  (int)
 *      Column in which row elements are currently being exchanged.
 *  Row1Ptr  (ElementPtr)
 *      Pointer to an element in Row1.
 *  Row2Ptr  (ElementPtr)
 *      Pointer to an element in Row2.
 *  Element1  (ElementPtr)
 *      Pointer to the element in Row1 to be exchanged.
 *  Element2  (ElementPtr)
 *      Pointer to the element in Row2 to be exchanged.
 */

spcRowExchange( Matrix, Row1, Row2 )

MatrixPtr Matrix;
int  Row1, Row2;
{
register  ElementPtr  Row1Ptr, Row2Ptr;
int  Column;
ElementPtr  Element1, Element2;

/* Begin `spcRowExchange'. */
    if (Row1 > Row2)  SWAP(int, Row1, Row2);

    Row1Ptr = Matrix->FirstInRow[Row1];
    Row2Ptr = Matrix->FirstInRow[Row2];
    while (Row1Ptr != NULL OR Row2Ptr != NULL)
    {
/* Exchange elements in rows while traveling from left to right. */
        if (Row1Ptr == NULL)
        {   Column = Row2Ptr->Col;
            Element1 = NULL;
            Element2 = Row2Ptr;
            Row2Ptr = Row2Ptr->NextInRow;
        }
        else if (Row2Ptr == NULL)
        {   Column = Row1Ptr->Col;
            Element1 = Row1Ptr;
            Element2 = NULL;
            Row1Ptr = Row1Ptr->NextInRow;
        }
        else if (Row1Ptr->Col < Row2Ptr->Col)
        {   Column = Row1Ptr->Col;
            Element1 = Row1Ptr;
            Element2 = NULL;
            Row1Ptr = Row1Ptr->NextInRow;
        }
        else if (Row1Ptr->Col > Row2Ptr->Col)
        {   Column = Row2Ptr->Col;
            Element1 = NULL;
            Element2 = Row2Ptr;
            Row2Ptr = Row2Ptr->NextInRow;
        }
        else   /* Row1Ptr->Col == Row2Ptr->Col */
        {   Column = Row1Ptr->Col;
            Element1 = Row1Ptr;
            Element2 = Row2Ptr;
            Row1Ptr = Row1Ptr->NextInRow;
            Row2Ptr = Row2Ptr->NextInRow;
        }

        ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column);
    }  /* end of while(Row1Ptr != NULL OR Row2Ptr != NULL) */

    if (Matrix->InternalVectorsAllocated)
        SWAP( int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]);
    SWAP( ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]);
    SWAP( int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]);
#if TRANSLATE
    Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row1] ] = Row1;
    Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row2] ] = Row2;
#endif

    return;
}









/*
 *  EXCHANGE COLUMNS
 *
 *  Performs all required operations to exchange two columns. Those operations
 *  include: swap FirstInCol pointers, fixing up the NextInRow pointers,
 *  swapping column indexes in MatrixElements, and swapping Markowitz 
 *  column counts.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  Col1  <input>  (int)
 *      Column index of one of the columns, becomes the smallest index.
 *  Col2  <input>  (int)
 *      Column index of the other column, becomes the largest index
 *
 *  Local variables:
 *  Row  (int)
 *      Row in which column elements are currently being exchanged.
 *  Col1Ptr  (ElementPtr)
 *      Pointer to an element in Col1.
 *  Col2Ptr  (ElementPtr)
 *      Pointer to an element in Col2.
 *  Element1  (ElementPtr)
 *      Pointer to the element in Col1 to be exchanged.
 *  Element2  (ElementPtr)
 *      Pointer to the element in Col2 to be exchanged.
 */

spcColExchange( Matrix, Col1, Col2 )

MatrixPtr Matrix;
int  Col1, Col2;
{
register  ElementPtr  Col1Ptr, Col2Ptr;
int  Row;
ElementPtr  Element1, Element2;

/* Begin `spcColExchange'. */
    if (Col1 > Col2)  SWAP(int, Col1, Col2);

    Col1Ptr = Matrix->FirstInCol[Col1];
    Col2Ptr = Matrix->FirstInCol[Col2];
    while (Col1Ptr != NULL OR Col2Ptr != NULL)
    {
/* Exchange elements in rows while traveling from top to bottom. */
        if (Col1Ptr == NULL)
        {   Row = Col2Ptr->Row;
            Element1 = NULL;
            Element2 = Col2Ptr;
            Col2Ptr = Col2Ptr->NextInCol;
        }
        else if (Col2Ptr == NULL)
        {   Row = Col1Ptr->Row;
            Element1 = Col1Ptr;
            Element2 = NULL;
            Col1Ptr = Col1Ptr->NextInCol;
        }
        else if (Col1Ptr->Row < Col2Ptr->Row)
        {   Row = Col1Ptr->Row;
            Element1 = Col1Ptr;
            Element2 = NULL;
            Col1Ptr = Col1Ptr->NextInCol;
        }
        else if (Col1Ptr->Row > Col2Ptr->Row)
        {   Row = Col2Ptr->Row;
            Element1 = NULL;
            Element2 = Col2Ptr;
            Col2Ptr = Col2Ptr->NextInCol;
        }
        else   /* Col1Ptr->Row == Col2Ptr->Row */
        {   Row = Col1Ptr->Row;
            Element1 = Col1Ptr;
            Element2 = Col2Ptr;
            Col1Ptr = Col1Ptr->NextInCol;
            Col2Ptr = Col2Ptr->NextInCol;
        }

        ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row);
    }  /* end of while(Col1Ptr != NULL OR Col2Ptr != NULL) */

    if (Matrix->InternalVectorsAllocated)
        SWAP( int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]);
    SWAP( ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
    SWAP( int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
#if TRANSLATE
    Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col1] ] = Col1;
    Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col2] ] = Col2;
#endif

    return;
}







/*
 *  EXCHANGE TWO ELEMENTS IN A COLUMN
 *
 *  Performs all required operations to exchange two elements in a column.
 *  Those operations are: restring NextInCol pointers and swapping row indexes
 *  in the MatrixElements.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  Row1  <input>  (int)
 *      Row of top element to be exchanged.
 *  Element1  <input>  (ElementPtr)
 *      Pointer to top element to be exchanged.
 *  Row2  <input>  (int)
 *      Row of bottom element to be exchanged.
 *  Element2  <input>  (ElementPtr)
 *      Pointer to bottom element to be exchanged.
 *  Column <input>  (int)
 *      Column that exchange is to take place in.
 *
 *  >>> Local variables:
 *  ElementAboveRow1  (ElementPtr *)
 *      Location of pointer which points to the element above Element1. This
 *      pointer is modified so that it points to correct element on exit.
 *  ElementAboveRow2  (ElementPtr *)
 *      Location of pointer which points to the element above Element2. This
 *      pointer is modified so that it points to correct element on exit.
 *  ElementBelowRow1  (ElementPtr)
 *      Pointer to element below Element1.
 *  ElementBelowRow2  (ElementPtr)
 *      Pointer to element below Element2.
 *  pElement  (ElementPtr)
 *      Pointer used to traverse the column.
 */

static
ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column )

MatrixPtr Matrix;
register  ElementPtr  Element1, Element2;
int  Row1, Row2, Column;
{
ElementPtr  *ElementAboveRow1, *ElementAboveRow2;
ElementPtr  ElementBelowRow1, ElementBelowRow2;
register  ElementPtr  pElement;

/* Begin `ExchangeColElements'. */
/* Search to find the ElementAboveRow1. */
    ElementAboveRow1 = &(Matrix->FirstInCol[Column]);
    pElement = *ElementAboveRow1;
    while (pElement->Row < Row1)
    {   ElementAboveRow1 = &(pElement->NextInCol);
        pElement = *ElementAboveRow1;
    }
    if (Element1 != NULL)
    {   ElementBelowRow1 = Element1->NextInCol;
        if (Element2 == NULL)
        {
/* Element2 does not exist, move Element1 down to Row2. */
            if ( ElementBelowRow1 != NULL AND ElementBelowRow1->Row < Row2 )
            {
/* Element1 must be removed from linked list and moved. */
                *ElementAboveRow1 = ElementBelowRow1;

/* Search column for Row2. */
                pElement = ElementBelowRow1;
                do
                {   ElementAboveRow2 = &(pElement->NextInCol);
                    pElement = *ElementAboveRow2;
                }   while (pElement != NULL AND pElement->Row < Row2);

/* Place Element1 in Row2. */
                *ElementAboveRow2 = Element1;
                Element1->NextInCol = pElement;
                *ElementAboveRow1 =ElementBelowRow1;
            }
            Element1->Row = Row2;
        }
        else
        {
/* Element2 does exist, and the two elements must be exchanged. */
            if ( ElementBelowRow1->Row == Row2)
            {
/* Element2 is just below Element1, exchange them. */
                Element1->NextInCol = Element2->NextInCol;
                Element2->NextInCol = Element1;
                *ElementAboveRow1 = Element2;
            }
            else
            {
/* Element2 is not just below Element1 and must be searched for. */
                pElement = ElementBelowRow1;
                do
                {   ElementAboveRow2 = &(pElement->NextInCol);
                    pElement = *ElementAboveRow2;
                }   while (pElement->Row < Row2);

                ElementBelowRow2 = Element2->NextInCol;

/* Switch Element1 and Element2. */
                *ElementAboveRow1 = Element2;
                Element2->NextInCol = ElementBelowRow1;
                *ElementAboveRow2 = Element1;
                Element1->NextInCol = ElementBelowRow2;
            }
            Element1->Row = Row2;
            Element2->Row = Row1;
        }
    }
    else
    {
/* Element1 does not exist. */
        ElementBelowRow1 = pElement;

/* Find Element2. */
        if (ElementBelowRow1->Row != Row2)
        {   do
            {   ElementAboveRow2 = &(pElement->NextInCol);
                pElement = *ElementAboveRow2;
            }   while (pElement->Row < Row2);

        ElementBelowRow2 = Element2->NextInCol;

/* Move Element2 to Row1. */
            *ElementAboveRow2 = Element2->NextInCol;
            *ElementAboveRow1 = Element2;
            Element2->NextInCol = ElementBelowRow1;
        }
        Element2->Row = Row1;
    }
    return;
}







/*
 *  EXCHANGE TWO ELEMENTS IN A ROW
 *
 *  Performs all required operations to exchange two elements in a row.
 *  Those operations are: restring NextInRow pointers and swapping column
 *  indexes in the MatrixElements.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  Col1  <input>  (int)
 *      Col of left-most element to be exchanged.
 *  Element1  <input>  (ElementPtr)
 *      Pointer to left-most element to be exchanged.
 *  Col2  <input>  (int)
 *      Col of right-most element to be exchanged.
 *  Element2  <input>  (ElementPtr)
 *      Pointer to right-most element to be exchanged.
 *  Row <input>  (int)
 *      Row that exchange is to take place in.
 *
 *  >>> Local variables:
 *  ElementLeftOfCol1  (ElementPtr *)
 *      Location of pointer which points to the element to the left of
 *      Element1. This pointer is modified so that it points to correct
 *      element on exit.
 *  ElementLeftOfCol2  (ElementPtr *)
 *      Location of pointer which points to the element to the left of
 *      Element2. This pointer is modified so that it points to correct
 *      element on exit.
 *  ElementRightOfCol1  (ElementPtr)
 *      Pointer to element right of Element1.
 *  ElementRightOfCol2  (ElementPtr)
 *      Pointer to element right of Element2.
 *  pElement  (ElementPtr)
 *      Pointer used to traverse the row.
 */

static
ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row )

MatrixPtr Matrix;
int  Col1, Col2, Row;
register ElementPtr  Element1, Element2;
{
ElementPtr  *ElementLeftOfCol1, *ElementLeftOfCol2;
ElementPtr  ElementRightOfCol1, ElementRightOfCol2;
register   ElementPtr  pElement;

/* Begin `ExchangeRowElements'. */
/* Search to find the ElementLeftOfCol1. */
    ElementLeftOfCol1 = &(Matrix->FirstInRow[Row]);
    pElement = *ElementLeftOfCol1;
    while (pElement->Col < Col1)
    {   ElementLeftOfCol1 = &(pElement->NextInRow);
        pElement = *ElementLeftOfCol1;
    }
    if (Element1 != NULL)
    {   ElementRightOfCol1 = Element1->NextInRow;
        if (Element2 == NULL)
        {
/* Element2 does not exist, move Element1 to right to Col2. */
            if ( ElementRightOfCol1 != NULL AND ElementRightOfCol1->Col < Col2 )
            {
/* Element1 must be removed from linked list and moved. */
                *ElementLeftOfCol1 = ElementRightOfCol1;

/* Search Row for Col2. */
                pElement = ElementRightOfCol1;
                do
                {   ElementLeftOfCol2 = &(pElement->NextInRow);
                    pElement = *ElementLeftOfCol2;
                }   while (pElement != NULL AND pElement->Col < Col2);

/* Place Element1 in Col2. */
                *ElementLeftOfCol2 = Element1;
                Element1->NextInRow = pElement;
                *ElementLeftOfCol1 =ElementRightOfCol1;
            }
            Element1->Col = Col2;
        }
        else
        {
/* Element2 does exist, and the two elements must be exchanged. */
            if ( ElementRightOfCol1->Col == Col2)
            {
/* Element2 is just right of Element1, exchange them. */
                Element1->NextInRow = Element2->NextInRow;
                Element2->NextInRow = Element1;
                *ElementLeftOfCol1 = Element2;
            }
            else
            {
/* Element2 is not just right of Element1 and must be searched for. */
                pElement = ElementRightOfCol1;
                do
                {   ElementLeftOfCol2 = &(pElement->NextInRow);
                    pElement = *ElementLeftOfCol2;
                }   while (pElement->Col < Col2);

                ElementRightOfCol2 = Element2->NextInRow;

/* Switch Element1 and Element2. */
                *ElementLeftOfCol1 = Element2;
                Element2->NextInRow = ElementRightOfCol1;
                *ElementLeftOfCol2 = Element1;
                Element1->NextInRow = ElementRightOfCol2;
            }
            Element1->Col = Col2;
            Element2->Col = Col1;
        }
    }
    else
    {
/* Element1 does not exist. */
        ElementRightOfCol1 = pElement;

/* Find Element2. */
        if (ElementRightOfCol1->Col != Col2)
        {   do
            {   ElementLeftOfCol2 = &(pElement->NextInRow);
                pElement = *ElementLeftOfCol2;
            }   while (pElement->Col < Col2);

            ElementRightOfCol2 = Element2->NextInRow;

/* Move Element2 to Col1. */
            *ElementLeftOfCol2 = Element2->NextInRow;
            *ElementLeftOfCol1 = Element2;
            Element2->NextInRow = ElementRightOfCol1;
        }
        Element2->Col = Col1;
    }
    return;
}











/*
 *  PERFORM ROW AND COLUMN ELIMINATION ON REAL MATRIX
 *
 *  Eliminates a single row and column of the matrix and leaves single row of
 *  the upper triangular matrix and a single column of the lower triangular
 *  matrix in its wake.  Uses Gauss's method.
 *
 *  >>> Argument:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  pPivot  <input>  (ElementPtr)
 *      Pointer to the current pivot.
 *
 *  >>> Local variables:
 *  pLower  (ElementPtr)
 *      Points to matrix element in lower triangular column.
 *  pSub (ElementPtr)
 *      Points to elements in the reduced submatrix.
 *  Row  (int)
 *      Row index.
 *  pUpper  (ElementPtr)
 *      Points to matrix element in upper triangular row.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static
RealRowColElimination( Matrix, pPivot )

MatrixPtr Matrix;
register  ElementPtr  pPivot;
{
#if REAL
register  ElementPtr  pSub;
register  int  Row;
register  ElementPtr  pLower, pUpper;
extern ElementPtr  CreateFillin();

/* Begin `RealRowColElimination'. */

/* Test for zero pivot. */
    if (ABS(pPivot->Real) == 0.0)
    {   (void)MatrixIsSingular( Matrix, pPivot->Row );
        return;
    }
    pPivot->Real = 1.0 / pPivot->Real;

    pUpper = pPivot->NextInRow;
    while (pUpper != NULL)
    {
/* Calculate upper triangular element. */
        pUpper->Real *= pPivot->Real;

        pSub = pUpper->NextInCol;
        pLower = pPivot->NextInCol;
        while (pLower != NULL)
        {   Row = pLower->Row;

/* Find element in row that lines up with current lower triangular element. */
            while (pSub != NULL AND pSub->Row < Row)
                pSub = pSub->NextInCol;

/* Test to see if desired element was not found, if not, create fill-in. */
            if (pSub == NULL OR pSub->Row > Row)
            {   pSub = CreateFillin( Matrix, Row, pUpper->Col );
                if (pSub == NULL)
                {   Matrix->Error = spNO_MEMORY;
                    return;
                }
            }
            pSub->Real -= pUpper->Real * pLower->Real;
            pSub = pSub->NextInCol;
            pLower = pLower->NextInCol;
        }
        pUpper = pUpper->NextInRow;
    }
    return;
#endif /* REAL */
}









/*
 *  PERFORM ROW AND COLUMN ELIMINATION ON COMPLEX MATRIX
 *
 *  Eliminates a single row and column of the matrix and leaves single row of
 *  the upper triangular matrix and a single column of the lower triangular
 *  matrix in its wake.  Uses Gauss's method.
 *
 *  >>> Argument:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  pPivot  <input>  (ElementPtr)
 *      Pointer to the current pivot.
 *
 *  >>> Local variables:
 *  pLower  (ElementPtr)
 *      Points to matrix element in lower triangular column.
 *  pSub (ElementPtr)
 *      Points to elements in the reduced submatrix.
 *  Row  (int)
 *      Row index.
 *  pUpper  (ElementPtr)
 *      Points to matrix element in upper triangular row.
 *
 *  Possible errors:
 *  spNO_MEMORY
 */

static
ComplexRowColElimination( Matrix, pPivot )

MatrixPtr Matrix;
register  ElementPtr  pPivot;
{
#if spCOMPLEX
register  ElementPtr  pSub;
register  int  Row;
register  ElementPtr  pLower, pUpper;
ElementPtr  CreateFillin();

/* Begin `ComplexRowColElimination'. */

/* Test for zero pivot. */
    if (ELEMENT_MAG(pPivot) == 0.0)
    {   (void)MatrixIsSingular( Matrix, pPivot->Row );
        return;
    }
    CMPLX_RECIPROCAL(*pPivot, *pPivot);

    pUpper = pPivot->NextInRow;
    while (pUpper != NULL)
    {
/* Calculate upper triangular element. */
/* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */
        CMPLX_MULT_ASSIGN(*pUpper, *pPivot);

        pSub = pUpper->NextInCol;
        pLower = pPivot->NextInCol;
        while (pLower != NULL)
        {   Row = pLower->Row;

/* Find element in row that lines up with current lower triangular element. */
            while (pSub != NULL AND pSub->Row < Row)
                pSub = pSub->NextInCol;

/* Test to see if desired element was not found, if not, create fill-in. */
            if (pSub == NULL OR pSub->Row > Row)
            {   pSub = CreateFillin( Matrix, Row, pUpper->Col );
                if (pSub == NULL)
                {   Matrix->Error = spNO_MEMORY;
                    return;
                }
            }

/* Cmplx expr: pElement -= *pUpper * pLower. */
            CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower);
            pSub = pSub->NextInCol;
            pLower = pLower->NextInCol;
        }
        pUpper = pUpper->NextInRow;
    }
    return;
#endif /* spCOMPLEX */
}





/*
 *  UPDATE MARKOWITZ NUMBERS
 *
 *  Updates the Markowitz numbers after a row and column have been eliminated.
 *  Also updates singleton count.
 *
 *  >>> Argument:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  pPivot  <input>  (ElementPtr)
 *      Pointer to the current pivot.
 *
 *  >>> Local variables:
 *  Row  (int)
 *      Row index.
 *  Col  (int)
 *      Column index.
 *  ColPtr  (ElementPtr)
 *      Points to matrix element in upper triangular column.
 *  RowPtr  (ElementPtr)
 *      Points to matrix element in lower triangular row.
 */

static
UpdateMarkowitzNumbers( Matrix, pPivot )

MatrixPtr Matrix;
ElementPtr  pPivot;
{
register  int  Row, Col;
register  ElementPtr  ColPtr, RowPtr;
register  int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol;
double Product;

/* Begin `UpdateMarkowitzNumbers'. */

/* Update Markowitz numbers. */
    for (ColPtr = pPivot->NextInCol; ColPtr != NULL; ColPtr = ColPtr->NextInCol)
    {   Row = ColPtr->Row;
        --MarkoRow[Row];

/* Form Markowitz product while being cautious of overflows. */
        if ((MarkoRow[Row] > LARGEST_SHORT_INTEGER AND MarkoCol[Row] != 0) OR
            (MarkoCol[Row] > LARGEST_SHORT_INTEGER AND MarkoRow[Row] != 0))
        {   Product = MarkoCol[Row] * MarkoRow[Row];
            if (Product >= LARGEST_LONG_INTEGER)
                Matrix->MarkowitzProd[Row] = LARGEST_LONG_INTEGER;
            else
                Matrix->MarkowitzProd[Row] = Product;
        }
        else Matrix->MarkowitzProd[Row] = MarkoRow[Row] * MarkoCol[Row];
        if (MarkoRow[Row] == 0)
            Matrix->Singletons++;
    }

    for (RowPtr = pPivot->NextInRow; RowPtr != NULL; RowPtr = RowPtr->NextInRow)
    {   Col = RowPtr->Col;
        --MarkoCol[Col];

/* Form Markowitz product while being cautious of overflows. */
        if ((MarkoRow[Col] > LARGEST_SHORT_INTEGER AND MarkoCol[Col] != 0) OR
            (MarkoCol[Col] > LARGEST_SHORT_INTEGER AND MarkoRow[Col] != 0))
        {   Product = MarkoCol[Col] * MarkoRow[Col];
            if (Product >= LARGEST_LONG_INTEGER)
                Matrix->MarkowitzProd[Col] = LARGEST_LONG_INTEGER;
            else
                Matrix->MarkowitzProd[Col] = Product;
        }
        else Matrix->MarkowitzProd[Col] = MarkoRow[Col] * MarkoCol[Col];
        if ((MarkoCol[Col] == 0) AND (MarkoRow[Col] != 0))
            Matrix->Singletons++;
    }
    return;
}








/*
 *  CREATE FILL-IN
 *
 *  This routine is used to create fill-ins and splice them into the
 *  matrix.
 *
 *  >>> Returns:
 *  Pointer to fill-in.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to the matrix.
 *  Col  <input>  (int)
 *      Column index for element.
 *  Row  <input>  (int)
 *      Row index for element.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      Pointer to an element in the matrix.
 *  ppElementAbove  (ElementPtr *)
 *      This contains the address of the pointer to the element just above the
 *      one being created. It is used to speed the search and it is updated
 *      with address of the created element.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static ElementPtr
CreateFillin( Matrix, Row, Col )

MatrixPtr Matrix;
register int  Row;
int  Col;
{
register  ElementPtr  pElement, *ppElementAbove;
ElementPtr  spcCreateElement();

/* Begin `CreateFillin'. */

/* Find Element above fill-in. */
    ppElementAbove = &Matrix->FirstInCol[Col];
    pElement = *ppElementAbove;
    while (pElement != NULL)
    {   if (pElement->Row < Row)
        {   ppElementAbove = &pElement->NextInCol;
            pElement = *ppElementAbove;
        }
        else break; /* while loop */
    }

/* End of search, create the element. */
    pElement = spcCreateElement( Matrix, Row, Col, ppElementAbove, YES );

/* Update Markowitz counts and products. */
    Matrix->MarkowitzProd[Row] = ++Matrix->MarkowitzRow[Row] *
                                   Matrix->MarkowitzCol[Row];
    if ((Matrix->MarkowitzRow[Row] == 1) AND (Matrix->MarkowitzCol[Row] != 0))
        Matrix->Singletons--;
    Matrix->MarkowitzProd[Col] = ++Matrix->MarkowitzCol[Col] *
                                   Matrix->MarkowitzRow[Col];
    if ((Matrix->MarkowitzRow[Col] != 0) AND (Matrix->MarkowitzCol[Col] == 1))
        Matrix->Singletons--;

    return pElement;
}








/*
 *  ZERO PIVOT ENCOUNTERED
 *
 *  This routine is called when a singular matrix is found.  It then
 *  records the current row and column and exits.
 *
 *  >>> Returned:
 *  The error code spSINGULAR or spZERO_DIAG is returned.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *  Step  <input>  (int)
 *      Index of diagonal that is zero.
 */

static int
MatrixIsSingular( Matrix, Step )

MatrixPtr  Matrix;
int  Step;
{
/* Begin `MatrixIsSingular'. */

    Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
    Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
    return (Matrix->Error = spSINGULAR);
}


static int
ZeroPivot( Matrix, Step )

MatrixPtr  Matrix;
int  Step;
{
/* Begin `ZeroPivot'. */

    Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
    Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
    return (Matrix->Error = spZERO_DIAG);
}






#if ANNOTATE == FULL

/*
 *
 *  WRITE STATUS
 *
 *  Write a summary of important variables to standard output.
 */

static
WriteStatus( Matrix, Step )

MatrixPtr Matrix;
int Step;
{
int  I;

/* Begin `WriteStatus'. */

    printf("Step = %1d   ", Step);
    printf("Pivot found at %1d,%1d using ", Matrix->PivotsOriginalRow,
                                            Matrix->PivotsOriginalCol);
    switch(Matrix->PivotSelectionMethod)
    {   case 's': printf("SearchForSingleton\n");  break;
        case 'q': printf("QuicklySearchDiagonal\n");  break;
        case 'd': printf("SearchDiagonal\n");  break;
        case 'e': printf("SearchEntireMatrix\n");  break;
    }

    printf("MarkowitzRow     = ");
    for (I = 1; I <= Matrix->Size; I++)
        printf("%2d  ", Matrix->MarkowitzRow[I]);
    printf("\n");

    printf("MarkowitzCol     = ");
    for (I = 1; I <= Matrix->Size; I++)
        printf("%2d  ", Matrix->MarkowitzCol[I]);
    printf("\n");

    printf("MarkowitzProduct = ");
    for (I = 1; I <= Matrix->Size; I++)
        printf("%2d  ", Matrix->MarkowitzProd[I]);
    printf("\n");

    printf("Singletons = %2d\n", Matrix->Singletons);

    printf("IntToExtRowMap     = ");
    for (I = 1; I <= Matrix->Size; I++)
        printf("%2d  ", Matrix->IntToExtRowMap[I]);
    printf("\n");

    printf("IntToExtColMap     = ");
    for (I = 1; I <= Matrix->Size; I++)
        printf("%2d  ", Matrix->IntToExtColMap[I]);
    printf("\n");

    printf("ExtToIntRowMap     = ");
    for (I = 1; I <= Matrix->ExtSize; I++)
        printf("%2d  ", Matrix->ExtToIntRowMap[I]);
    printf("\n");

    printf("ExtToIntColMap     = ");
    for (I = 1; I <= Matrix->ExtSize; I++)
        printf("%2d  ", Matrix->ExtToIntColMap[I]);
    printf("\n\n");

/*  spPrint((char *)Matrix, NO, YES);  */

    return;

}
#endif /* ANNOTATE == FULL */
@EOF

chmod 644 spFactor.c

echo x - spMatrix.h
cat >spMatrix.h <<'@EOF'
/*
 *  EXPORTS for sparse matrix routines.
 *
 *  Author:                     Advising professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      UC Berkeley
 *
 *  This file contains definitions that are useful to the calling
 *  program.  In particular, this file contains error keyword
 *  definitions, some macro functions that are used to quickly enter
 *  data into the matrix and the type definition of a data structure
 *  that acts as a template for entering admittances into the matrix.
 *  Also included is the type definitions for the various functions
 *  available to the user.
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 *
 *  $Date: 88/06/18 11:14:02 $
 *  $Revision: 1.2 $
 */




#ifndef  spOKAY

/*
 *  IMPORTS
 *
 *  >>> Import descriptions:
 *  spConfig.h
 *      Macros that customize the sparse matrix routines.
 */

#include "spConfig.h"






/*
 *  ERROR KEYWORDS
 *
 *  The actual numbers used in the error codes are not sacred, they can be
 *  changed under the condition that the codes for the nonfatal errors are
 *  less than the code for spFATAL and similarly the codes for the fatal
 *  errors are greater than that for spFATAL.
 *
 *  >>> Error descriptions:
 *  spOKAY
 *      No error has occurred.
 *  spSMALL_PIVOT
 *      When reordering the matrix, no element was found which satisfies the
 *      absolute threshold criteria.  The largest element in the matrix was
 *      chosen as pivot.  Non-fatal.
 *  spZERO_DIAG
 *      Fatal error.  A zero was encountered on the diagonal the matrix.  This
 *      does not necessarily imply that the matrix is singular.  When this
 *      error occurs, the matrix should be reconstructed and factored using
 *      spOrderAndFactor().
 *  spSINGULAR
 *      Fatal error.  Matrix is singular, so no unique solution exists.
 *  spNO_MEMORY
 *      Fatal error.  Indicates that not enough memory is available to handle
 *      the matrix.
 *  spPANIC
 *      Fatal error indicating that the routines are not prepared to
 *      handle the matrix that has been requested.  This may occur when
 *      the matrix is specified to be real and the routines are not
 *      compiled for real matrices, or when the matrix is specified to
 *      be complex and the routines are not compiled to handle complex
 *      matrices.
 *  spFATAL
 *      Not an error flag, but rather the dividing line between fatal errors
 *      and warnings.
 */

/* Begin error macros. */
#define  spOKAY                 0
#define  spSMALL_PIVOT          1
#define  spZERO_DIAG            2
#define  spSINGULAR             3
#define  spNO_MEMORY            4
#define  spPANIC                5

#define  spFATAL                2







/*
 *  KEYWORD DEFINITIONS
 *
 *  Here we define what precision arithmetic Sparse will use.  Double
 *  precision is suggested as being most appropriate for circuit
 *  simulation and for C.  However, it is possible to change spREAL
 *  to a float for single precision arithmetic.  Note that in C, single
 *  precision arithmetic is often slower than double precision.  Sparse
 *  internally refers to spREALs as RealNumbers.
 *
 *  Some C compilers, notably the old VMS compiler, do not handle the keyword
 *  "void" correctly.  If this is true for your compiler, remove the
 *  comment delimiters from the redefinition of void to int below.
 */

#define  spREAL double
/* #define  void    int   */





/*
 *  PARTITION TYPES
 *
 *  When factoring a previously ordered matrix using spFactor(), Sparse
 *  operates on a row-at-a-time basis.  For speed, on each step, the row
 *  being updated is copied into a full vector and the operations are
 *  performed on that vector.  This can be done one of two ways, either
 *  using direct addressing or indirect addressing.  Direct addressing
 *  is fastest when the matrix is relatively dense and indirect addressing
 *  is quite sparse.  The user can select which partitioning mode is used.
 *  The following keywords are passed to spPartition() and indicate that
 *  Sparse should use only direct addressing, only indirect addressing, or
 *  that it should choose the best mode on a row-by-row basis.  The time
 *  required to choose a partition is of the same order of the cost to factor
 *  the matrix.
 *
 *  If you plan to factor a large number of matrices with the same structure,
 *  it is best to let Sparse choose the partition.  Otherwise, you should
 *  choose the partition based on the predicted density of the matrix.
 */

/* Begin partition keywords. */

#define spDEFAULT_PARTITION     0
#define spDIRECT_PARTITION      1
#define spINDIRECT_PARTITION    2
#define spAUTO_PARTITION        3





/*
 *  MACRO FUNCTION DEFINITIONS
 *
 *  >>> Macro descriptions:
 *  spADD_REAL_ELEMENT
 *      Macro function that adds data to a real element in the matrix by a
 *      pointer.
 *  spADD_IMAG_ELEMENT
 *      Macro function that adds data to a imaginary element in the matrix by
 *      a pointer.
 *  spADD_COMPLEX_ELEMENT
 *      Macro function that adds data to a complex element in the matrix by a
 *      pointer.
 *  spADD_REAL_QUAD
 *      Macro function that adds data to each of the four real matrix elements
 *      specified by the given template.
 *  spADD_IMAG_QUAD
 *      Macro function that adds data to each of the four imaginary matrix
 *      elements specified by the given template.
 *  spADD_COMPLEX_QUAD
 *      Macro function that adds data to each of the four complex matrix
 *      elements specified by the given template.
 */

/* Begin Macros. */
#define  spADD_REAL_ELEMENT(element,real)       *(element) += real

#define  spADD_IMAG_ELEMENT(element,imag)       *(element+1) += imag

#define  spADD_COMPLEX_ELEMENT(element,real,imag)       \
{   *(element) += real;                                 \
    *(element+1) += imag;                               \
}

#define  spADD_REAL_QUAD(template,real)         \
{   *((template).Element1) += real;             \
    *((template).Element2) += real;             \
    *((template).Element3Negated) -= real;      \
    *((template).Element4Negated) -= real;      \
}

#define  spADD_IMAG_QUAD(template,imag)         \
{   *((template).Element1+1) += imag;           \
    *((template).Element2+1) += imag;           \
    *((template).Element3Negated+1) -= imag;    \
    *((template).Element4Negated+1) -= imag;    \
}

#define  spADD_COMPLEX_QUAD(template,real,imag) \
{   *((template).Element1) += real;             \
    *((template).Element2) += real;             \
    *((template).Element3Negated) -= real;      \
    *((template).Element4Negated) -= real;      \
    *((template).Element1+1) += imag;           \
    *((template).Element2+1) += imag;           \
    *((template).Element3Negated+1) -= imag;    \
    *((template).Element4Negated+1) -= imag;    \
}







/*
 *   TYPE DEFINITION FOR COMPONENT TEMPLATE
 *
 *   This data structure is used to hold pointers to four related elements in
 *   matrix.  It is used in conjunction with the routines
 *       spGetAdmittance
 *       spGetQuad
 *       spGetOnes
 *   These routines stuff the structure which is later used by the spADD_QUAD
 *   macro functions above.  It is also possible for the user to collect four
 *   pointers returned by spGetElement and stuff them into the template.
 *   The spADD_QUAD routines stuff data into the matrix in locations specified
 *   by Element1 and Element2 without changing the data.  The data is negated
 *   before being placed in Element3 and Element4.
 */

/* Begin `spTemplate'. */
struct  spTemplate
{   spREAL    *Element1       ;
    spREAL    *Element2       ;
    spREAL    *Element3Negated;
    spREAL    *Element4Negated;
};





/*
 *   FUNCTION TYPE DEFINITIONS
 *
 *   The type of every user accessible function is declared here.
 */

/* Begin function declarations. */

#ifdef __STDC__

/* For compilers that understand function prototypes. */

extern  void     spClear( char* );
extern  spREAL   spCondition( char*, spREAL, int* );
extern  char    *spCreate( int, int, int* );
extern  void     spDeleteRowAndCol( char*, int, int );
extern  void     spDestroy( char* );
extern  int      spElementCount( char* );
extern  int      spError( char* );
extern  int      spFactor( char* );
extern  int      spFileMatrix( char*, char*, char*, int, int, int );
extern  int      spFileStats( char*, char*, char* );
extern  int      spFillinCount( char* );
extern  int      spGetAdmittance( char*, int, int, struct spTemplate* );
extern  spREAL  *spGetElement( char*, int, int );
extern  char    *spGetInitInfo( spREAL* );
extern  int      spGetOnes( char*, int, int, int, struct spTemplate* );
extern  int      spGetQuad( char*, int, int, int, int, struct spTemplate* );
extern  int      spGetSize( char*, int );
extern  int      spInitialize( char*, int (*)() );
extern  void     spInstallInitInfo( spREAL*, char* );
extern  spREAL   spLargestElement( char* );
extern  void     spMNA_Preorder( char* );
extern  spREAL   spNorm( char* );
extern  int      spOrderAndFactor( char*, spREAL[], spREAL, spREAL, int );
extern  void     spPartition( char*, int );
extern  void     spPrint( char*, int, int, int );
extern  spREAL   spPseudoCondition( char* );
extern  spREAL   spRoundoff( char*, spREAL );
extern  void     spScale( char*, spREAL[], spREAL[] );
extern  void     spSetComplex( char* );
extern  void     spSetReal( char* );
extern  void     spStripFills( char* );
extern  void     spWhereSingular( char*, int*, int* );

/* Functions with argument lists that are dependent on options. */

#if spCOMPLEX
extern  void     spDeterminant ( char*, int*, spREAL*, spREAL* );
#else /* NOT spCOMPLEX */
extern  void     spDeterminant ( char*, int*, spREAL* );
#endif /* NOT spCOMPLEX */
#if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS
extern  int      spFileVector( char*, char* , spREAL[], spREAL[]);
extern  void     spMultiply( char*, spREAL[], spREAL[], spREAL[], spREAL[] );
extern  void     spMultTransposed(char*,spREAL[],spREAL[],spREAL[],spREAL[]);
extern  void     spSolve( char*, spREAL[], spREAL[], spREAL[], spREAL[] );
extern  void     spSolveTransposed(char*,spREAL[],spREAL[],spREAL[],spREAL[]);
#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */
extern  int      spFileVector( char*, char* , spREAL[] );
extern  void     spMultiply( char*, spREAL[], spREAL[] );
extern  void     spMultTransposed( char*, spREAL[], spREAL[] );
extern  void     spSolve( char*, spREAL[], spREAL[] );
extern  void     spSolveTransposed( char*, spREAL[], spREAL[] );
#endif /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */

#else /* NOT defined(__STDC__) */

/* For compilers that do not understand function prototypes. */

extern  void     spClear();
extern  spREAL   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  spREAL  *spGetElement();
extern  char    *spGetInitInfo();
extern  int      spGetOnes();
extern  int      spGetQuad();
extern  int      spGetSize();
extern  int      spInitialize();
extern  void     spInstallInitInfo();
extern  spREAL   spLargestElement();
extern  void     spMNA_Preorder();
extern  void     spMultiply();
extern  void     spMultTransposed();
extern  spREAL   spNorm();
extern  int      spOrderAndFactor();
extern  void     spPartition();
extern  void     spPrint();
extern  spREAL   spPseudoCondition();
extern  spREAL   spRoundoff();
extern  void     spScale();
extern  void     spSetComplex();
extern  void     spSetReal();
extern  void     spSolve();
extern  void     spSolveTransposed();
extern  void     spStripFills();
extern  void     spWhereSingular();
#endif /* defined(__STDC__) */

#endif  /* spOKAY */
@EOF

chmod 644 spMatrix.h

echo x - spSolve.c
cat >spSolve.c <<'@EOF'
/*
 *  MATRIX SOLVE MODULE
 *
 *  Author:                     Advising professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      UC Berkeley
 *
 *  This file contains the forward and backward substitution routines for
 *  the sparse matrix routines.
 *
 *  >>> User accessible functions contained in this file:
 *  spSolve
 *  spSolveTransposed
 *
 *  >>> Other functions contained in this file:
 *  SolveComplexMatrix
 *  SolveComplexTransposedMatrix
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  by Kenneth S. Kundert and the University of California.
 *
 *  Permission to use, copy, modify, and distribute this software and
 *  its documentation for any purpose and without fee is hereby granted,
 *  provided that the copyright notices appear in all copies and
 *  supporting documentation and that the authors and the University of
 *  California are properly credited.  The authors and the University of
 *  California make no representations as to the suitability of this
 *  software for any purpose.  It is provided `as is', without express
 *  or implied warranty.
 */

#ifndef lint
static char copyright[] =
    "Sparse1.3: Copyright (c) 1985,86,87,88 by Kenneth S. Kundert";
static char RCSid[] =
    "@(#)$Header: spSolve.c,v 1.2 88/06/18 11:15:54 kundert Exp $";
#endif



/*
 *  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"





/*
 *  SOLVE MATRIX EQUATION
 *
 *  Performs forward elimination and back substitution to find the
 *  unknown vector from the RHS vector and factored matrix.  This
 *  routine assumes that the pivots are associated with the lower
 *  triangular (L) matrix and that the diagonal of the upper triangular
 *  (U) matrix consists of ones.  This routine arranges the computation
 *  in different way than is traditionally used in order to exploit the
 *  sparsity of the right-hand side.  See the reference in spRevision.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *  RHS  <input>  (RealVector)
 *      RHS is the input data array, the right hand side. This data is
 *      undisturbed and may be reused for other solves.
 *  Solution  <output>  (RealVector)
 *      Solution is the output data array. This routine is constructed such that
 *      RHS and Solution can be the same array.
 *  iRHS  <input>  (RealVector)
 *      iRHS is the imaginary portion of the input data array, the right
 *      hand side. This data is undisturbed and may be reused for other solves.
 *      This argument is only necessary if matrix is complex and if
 *      spSEPARATED_COMPLEX_VECTOR is set true.
 *  iSolution  <output>  (RealVector)
 *      iSolution is the imaginary portion of the output data array. This
 *      routine is constructed such that iRHS and iSolution can be
 *      the same array.  This argument is only necessary if matrix is complex
 *      and if spSEPARATED_COMPLEX_VECTOR is set true.
 *
 *  >>> Local variables:
 *  Intermediate  (RealVector)
 *      Temporary storage for use in forward elimination and backward
 *      substitution.  Commonly referred to as c, when the LU factorization
 *      equations are given as  Ax = b, Lc = b, Ux = c Local version of
 *      Matrix->Intermediate, which was created during the initial
 *      factorization in function CreateInternalVectors() in the matrix
 *      factorization module.
 *  pElement  (ElementPtr)
 *      Pointer used to address elements in both the lower and upper triangle
 *      matrices.
 *  pExtOrder  (int *)
 *      Pointer used to sequentially access each entry in IntToExtRowMap
 *      and IntToExtColMap arrays.  Used to quickly scramble and unscramble
 *      RHS and Solution to account for row and column interchanges.
 *  pPivot  (ElementPtr)
 *      Pointer that points to current pivot or diagonal element.
 *  Size  (int)
 *      Size of matrix. Made local to reduce indirection.
 *  Temp  (RealNumber)
 *      Temporary storage for entries in arrays.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

/*VARARGS3*/

void
spSolve( eMatrix, RHS, Solution IMAG_VECTORS )

char *eMatrix;
RealVector  RHS, Solution IMAG_VECTORS;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement;
register  RealVector  Intermediate;
register  RealNumber  Temp;
register  int  I, *pExtOrder, Size;
ElementPtr  pPivot;
void SolveComplexMatrix();

/* Begin `spSolve'. */
    ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) );

#if spCOMPLEX
    if (Matrix->Complex)
    {   SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS );
        return;
    }
#endif

#if REAL
    Intermediate = Matrix->Intermediate;
    Size = Matrix->Size;

/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
    --RHS;
    --Solution;
#endif

/* Initialize Intermediate vector. */
    pExtOrder = &Matrix->IntToExtRowMap[Size];
    for (I = Size; I > 0; I--)
        Intermediate[I] = RHS[*(pExtOrder--)];

/* Forward elimination. Solves Lc = b.*/
    for (I = 1; I <= Size; I++)
    {   
/* This step of the elimination is skipped if Temp equals zero. */
        if ((Temp = Intermediate[I]) != 0.0)
        {   pPivot = Matrix->Diag[I];
            Intermediate[I] = (Temp *= pPivot->Real);

            pElement = pPivot->NextInCol;
            while (pElement != NULL)
            {   Intermediate[pElement->Row] -= Temp * pElement->Real;
                pElement = pElement->NextInCol;
            }
        }
    }

/* Backward Substitution. Solves Ux = c.*/
    for (I = Size; I > 0; I--)
    {   Temp = Intermediate[I];
        pElement = Matrix->Diag[I]->NextInRow;
        while (pElement != NULL)
        {   Temp -= pElement->Real * Intermediate[pElement->Col];
            pElement = pElement->NextInRow;
        }
        Intermediate[I] = Temp;
    }

/* Unscramble Intermediate vector while placing data in to Solution vector. */
    pExtOrder = &Matrix->IntToExtColMap[Size];
    for (I = Size; I > 0; I--)
        Solution[*(pExtOrder--)] = Intermediate[I];

    return;
#endif /* REAL */
}











#if spCOMPLEX
/*
 *  SOLVE COMPLEX MATRIX EQUATION
 *
 *  Performs forward elimination and back substitution to find the
 *  unknown vector from the RHS vector and factored matrix.  This
 *  routine assumes that the pivots are associated with the lower
 *  triangular (L) matrix and that the diagonal of the upper triangular
 *  (U) matrix consists of ones.  This routine arranges the computation
 *  in different way than is traditionally used in order to exploit the
 *  sparsity of the right-hand side.  See the reference in spRevision.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *  RHS  <input>  (RealVector)
 *      RHS is the real portion of the input data array, the right hand
 *      side. This data is undisturbed and may be reused for other solves.
 *  Solution  <output>  (RealVector)
 *      Solution is the real portion of the output data array. This routine
 *      is constructed such that RHS and Solution can be the same
 *      array.
 *  iRHS  <input>  (RealVector)
 *      iRHS is the imaginary portion of the input data array, the right
 *      hand side. This data is undisturbed and may be reused for other solves.
 *      If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to
 *      supply this array.
 *  iSolution  <output>  (RealVector)
 *      iSolution is the imaginary portion of the output data array. This
 *      routine is constructed such that iRHS and iSolution can be
 *      the same array.  If spSEPARATED_COMPLEX_VECTOR is set false, there is no
 *      need to supply this array.
 *
 *  >>> Local variables:
 *  Intermediate  (ComplexVector)
 *      Temporary storage for use in forward elimination and backward
 *      substitution.  Commonly referred to as c, when the LU factorization
 *      equations are given as  Ax = b, Lc = b, Ux = c.
 *      Local version of Matrix->Intermediate, which was created during
 *      the initial factorization in function CreateInternalVectors() in the
 *      matrix factorization module.
 *  pElement  (ElementPtr)
 *      Pointer used to address elements in both the lower and upper triangle
 *      matrices.
 *  pExtOrder  (int *)
 *      Pointer used to sequentially access each entry in IntToExtRowMap
 *      and IntToExtColMap arrays.  Used to quickly scramble and unscramble
 *      RHS and Solution to account for row and column interchanges.
 *  pPivot  (ElementPtr)
 *      Pointer that points to current pivot or diagonal element.
 *  Size  (int)
 *      Size of matrix. Made local to reduce indirection.
 *  Temp  (ComplexNumber)
 *      Temporary storage for entries in arrays.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

static void
SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS )

MatrixPtr  Matrix;
RealVector  RHS, Solution IMAG_VECTORS;
{
register  ElementPtr  pElement;
register  ComplexVector  Intermediate;
register  int  I, *pExtOrder, Size;
ElementPtr  pPivot;
register ComplexVector  ExtVector;
ComplexNumber  Temp;

/* Begin `SolveComplexMatrix'. */

    Size = Matrix->Size;
    Intermediate = (ComplexVector)Matrix->Intermediate;

/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
    --RHS;      --iRHS;
    --Solution; --iSolution;
#else
    RHS -= 2; Solution -= 2;
#endif
#endif

/* Initialize Intermediate vector. */
    pExtOrder = &Matrix->IntToExtRowMap[Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Size; I > 0; I--)
    {   Intermediate[I].Real = RHS[*(pExtOrder)];
        Intermediate[I].Imag = iRHS[*(pExtOrder--)];
    }
#else
    ExtVector = (ComplexVector)RHS;
    for (I = Size; I > 0; I--)
        Intermediate[I] = ExtVector[*(pExtOrder--)];
#endif

/* Forward substitution. Solves Lc = b.*/
    for (I = 1; I <= Size; I++)
    {   Temp = Intermediate[I];

/* This step of the substitution is skipped if Temp equals zero. */
        if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0))
        {   pPivot = Matrix->Diag[I];
/* Cmplx expr: Temp *= (1.0 / Pivot). */
            CMPLX_MULT_ASSIGN(Temp, *pPivot);
            Intermediate[I] = Temp;
            pElement = pPivot->NextInCol;
            while (pElement != NULL)
            {
/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
                CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row],
                                       Temp, *pElement);
                pElement = pElement->NextInCol;
            }
        }
    }

/* Backward Substitution. Solves Ux = c.*/
    for (I = Size; I > 0; I--)
    {   Temp = Intermediate[I];
        pElement = Matrix->Diag[I]->NextInRow;

        while (pElement != NULL)
        {
/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */
            CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]);
            pElement = pElement->NextInRow;
        }
        Intermediate[I] = Temp;
    }

/* Unscramble Intermediate vector while placing data in to Solution vector. */
    pExtOrder = &Matrix->IntToExtColMap[Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Size; I > 0; I--)
    {   Solution[*(pExtOrder)] = Intermediate[I].Real;
        iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
    }
#else
    ExtVector = (ComplexVector)Solution;
    for (I = Size; I > 0; I--)
        ExtVector[*(pExtOrder--)] = Intermediate[I];
#endif

    return;
}
#endif /* spCOMPLEX */














#if TRANSPOSE
/*
 *  SOLVE TRANSPOSED MATRIX EQUATION
 *
 *  Performs forward elimination and back substitution to find the
 *  unknown vector from the RHS vector and transposed factored
 *  matrix. This routine is useful when performing sensitivity analysis
 *  on a circuit using the adjoint method.  This routine assumes that
 *  the pivots are associated with the untransposed lower triangular
 *  (L) matrix and that the diagonal of the untransposed upper
 *  triangular (U) matrix consists of ones.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *  RHS  <input>  (RealVector)
 *      RHS is the input data array, the right hand side. This data is
 *      undisturbed and may be reused for other solves.
 *  Solution  <output>  (RealVector)
 *      Solution is the output data array. This routine is constructed such that
 *      RHS and Solution can be the same array.
 *  iRHS  <input>  (RealVector)
 *      iRHS is the imaginary portion of the input data array, the right
 *      hand side. This data is undisturbed and may be reused for other solves.
 *      If spSEPARATED_COMPLEX_VECTOR is set false, or if matrix is real, there
 *      is no need to supply this array.
 *  iSolution  <output>  (RealVector)
 *      iSolution is the imaginary portion of the output data array. This
 *      routine is constructed such that iRHS and iSolution can be
 *      the same array.  If spSEPARATED_COMPLEX_VECTOR is set false, or if
 *      matrix is real, there is no need to supply this array.
 *
 *  >>> Local variables:
 *  Intermediate  (RealVector)
 *      Temporary storage for use in forward elimination and backward
 *      substitution.  Commonly referred to as c, when the LU factorization
 *      equations are given as  Ax = b, Lc = b, Ux = c.  Local version of
 *      Matrix->Intermediate, which was created during the initial
 *      factorization in function CreateInternalVectors() in the matrix
 *      factorization module.
 *  pElement  (ElementPtr)
 *      Pointer used to address elements in both the lower and upper triangle
 *      matrices.
 *  pExtOrder  (int *)
 *      Pointer used to sequentially access each entry in IntToExtRowMap
 *      and IntToExtRowMap arrays.  Used to quickly scramble and unscramble
 *      RHS and Solution to account for row and column interchanges.
 *  pPivot  (ElementPtr)
 *      Pointer that points to current pivot or diagonal element.
 *  Size  (int)
 *      Size of matrix. Made local to reduce indirection.
 *  Temp  (RealNumber)
 *      Temporary storage for entries in arrays.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

/*VARARGS3*/

void
spSolveTransposed( eMatrix, RHS, Solution IMAG_VECTORS )

char *eMatrix;
RealVector  RHS, Solution IMAG_VECTORS;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement;
register  RealVector  Intermediate;
register  int  I, *pExtOrder, Size;
ElementPtr  pPivot;
RealNumber  Temp;
void SolveComplexTransposedMatrix();

/* Begin `spSolveTransposed'. */
    ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) );

#if spCOMPLEX
    if (Matrix->Complex)
    {   SolveComplexTransposedMatrix( Matrix, RHS, Solution IMAG_VECTORS );
        return;
    }
#endif

#if REAL
    Size = Matrix->Size;
    Intermediate = Matrix->Intermediate;

/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
    --RHS;
    --Solution;
#endif

/* Initialize Intermediate vector. */
    pExtOrder = &Matrix->IntToExtColMap[Size];
    for (I = Size; I > 0; I--)
        Intermediate[I] = RHS[*(pExtOrder--)];

/* Forward elimination. */
    for (I = 1; I <= Size; I++)
    {   
/* This step of the elimination is skipped if Temp equals zero. */
        if ((Temp = Intermediate[I]) != 0.0)
        {   pElement = Matrix->Diag[I]->NextInRow;
            while (pElement != NULL)
            {   Intermediate[pElement->Col] -= Temp * pElement->Real;
                pElement = pElement->NextInRow;
            }

        }
    }

/* Backward Substitution. */
    for (I = Size; I > 0; I--)
    {   pPivot = Matrix->Diag[I];
        Temp = Intermediate[I];
        pElement = pPivot->NextInCol;
        while (pElement != NULL)
        {   Temp -= pElement->Real * Intermediate[pElement->Row];
            pElement = pElement->NextInCol;
        }
        Intermediate[I] = Temp * pPivot->Real;
    }

/* Unscramble Intermediate vector while placing data in to Solution vector. */
    pExtOrder = &Matrix->IntToExtRowMap[Size];
    for (I = Size; I > 0; I--)
        Solution[*(pExtOrder--)] = Intermediate[I];

    return;
#endif /* REAL */
}
#endif /* TRANSPOSE */










#if TRANSPOSE AND spCOMPLEX
/*
 *  SOLVE COMPLEX TRANSPOSED MATRIX EQUATION
 *
 *  Performs forward elimination and back substitution to find the
 *  unknown vector from the RHS vector and transposed factored
 *  matrix. This routine is useful when performing sensitivity analysis
 *  on a circuit using the adjoint method.  This routine assumes that
 *  the pivots are associated with the untransposed lower triangular
 *  (L) matrix and that the diagonal of the untransposed upper
 *  triangular (U) matrix consists of ones.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to matrix.
 *  RHS  <input>  (RealVector)
 *      RHS is the input data array, the right hand
 *      side. This data is undisturbed and may be reused for other solves.
 *      This vector is only the real portion if the matrix is complex and
 *      spSEPARATED_COMPLEX_VECTORS is set true.
 *  Solution  <output>  (RealVector)
 *      Solution is the real portion of the output data array. This routine
 *      is constructed such that RHS and Solution can be the same array.
 *      This vector is only the real portion if the matrix is complex and
 *      spSEPARATED_COMPLEX_VECTORS is set true.
 *  iRHS  <input>  (RealVector)
 *      iRHS is the imaginary portion of the input data array, the right
 *      hand side. This data is undisturbed and may be reused for other solves.
 *      If either spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set false, there
 *      is no need to supply this array.
 *  iSolution  <output>  (RealVector)
 *      iSolution is the imaginary portion of the output data array. This
 *      routine is constructed such that iRHS and iSolution can be
 *      the same array.  If spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set
 *      false, there is no need to supply this array.
 *
 *  >>> Local variables:
 *  Intermediate  (ComplexVector)
 *      Temporary storage for use in forward elimination and backward
 *      substitution.  Commonly referred to as c, when the LU factorization
 *      equations are given as  Ax = b, Lc = b, Ux = c.  Local version of
 *      Matrix->Intermediate, which was created during
 *      the initial factorization in function CreateInternalVectors() in the
 *      matrix factorization module.
 *  pElement  (ElementPtr)
 *      Pointer used to address elements in both the lower and upper triangle
 *      matrices.
 *  pExtOrder  (int *)
 *      Pointer used to sequentially access each entry in IntToExtRowMap
 *      and IntToExtColMap arrays.  Used to quickly scramble and unscramble
 *      RHS and Solution to account for row and column interchanges.
 *  pPivot  (ElementPtr)
 *      Pointer that points to current pivot or diagonal element.
 *  Size  (int)
 *      Size of matrix. Made local to reduce indirection.
 *  Temp  (ComplexNumber)
 *      Temporary storage for entries in arrays.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

static void
SolveComplexTransposedMatrix(Matrix, RHS, Solution IMAG_VECTORS )

MatrixPtr  Matrix;
RealVector  RHS, Solution IMAG_VECTORS;
{
register  ElementPtr  pElement;
register  ComplexVector  Intermediate;
register  int  I, *pExtOrder, Size;
register  ComplexVector  ExtVector;
ElementPtr  pPivot;
ComplexNumber  Temp;

/* Begin `SolveComplexTransposedMatrix'. */

    Size = Matrix->Size;
    Intermediate = (ComplexVector)Matrix->Intermediate;

/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
    --RHS;      --iRHS;
    --Solution; --iSolution;
#else
    RHS -= 2;   Solution -= 2;
#endif
#endif

/* Initialize Intermediate vector. */
    pExtOrder = &Matrix->IntToExtColMap[Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Size; I > 0; I--)
    {   Intermediate[I].Real = RHS[*(pExtOrder)];
        Intermediate[I].Imag = iRHS[*(pExtOrder--)];
    }
#else
    ExtVector = (ComplexVector)RHS;
    for (I = Size; I > 0; I--)
        Intermediate[I] = ExtVector[*(pExtOrder--)];
#endif

/* Forward elimination. */
    for (I = 1; I <= Size; I++)
    {   Temp = Intermediate[I];

/* This step of the elimination is skipped if Temp equals zero. */
        if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0))
        {   pElement = Matrix->Diag[I]->NextInRow;
            while (pElement != NULL)
            {
/* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */
                CMPLX_MULT_SUBT_ASSIGN( Intermediate[pElement->Col],
                                        Temp, *pElement);
                pElement = pElement->NextInRow;
            }
        }
    }

/* Backward Substitution. */
    for (I = Size; I > 0; I--)
    {   pPivot = Matrix->Diag[I];
        Temp = Intermediate[I];
        pElement = pPivot->NextInCol;

        while (pElement != NULL)
        {
/* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */
            CMPLX_MULT_SUBT_ASSIGN(Temp,Intermediate[pElement->Row],*pElement);

            pElement = pElement->NextInCol;
        }
/* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */
        CMPLX_MULT(Intermediate[I], Temp, *pPivot);
    }

/* Unscramble Intermediate vector while placing data in to Solution vector. */
    pExtOrder = &Matrix->IntToExtRowMap[Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Size; I > 0; I--)
    {   Solution[*(pExtOrder)] = Intermediate[I].Real;
        iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
    }
#else
    ExtVector = (ComplexVector)Solution;
    for (I = Size; I > 0; I--)
        ExtVector[*(pExtOrder--)] = Intermediate[I];
#endif

    return;
}
#endif /* TRANSPOSE AND spCOMPLEX */
@EOF

chmod 644 spSolve.c

exit 0


.