spfactor.c

ngspice又一个电子CAD仿真软件代码.功能更全

                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. */
            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);
}






/*
 *  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;
{
     ElementPtr  pElement;
     ElementPtr  pColumn;
     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. */
             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. */
             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);
}






/*
 *  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(void *eMatrix, int Mode)
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    ElementPtr  pElement, pColumn;
    int  Step, Size;
    int  *Nc, *No, *Nm;
    int *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++) {
            DoRealDirect[Step] = YES;
	    DoCmplxDirect[Step] = YES;
	}
        return;
    }
    else if (Mode == spINDIRECT_PARTITION)
    {
	for (Step = 1; Step <= Size; Step++)
	{
            DoRealDirect[Step] = NO;
	    DoCmplxDirect[Step] = NO;
	}
        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
        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
        /* On the hp350, it is never profitable to use direct for complex. */
        DoCmplxDirect[Step] = NO;
#undef generic
#endif

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

#ifdef generic
        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
#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
 */

void
spcCreateInternalVectors( Matrix )

MatrixPtr  Matrix;
{
    int  Size;

    /* Begin `spcCreateInternalVectors'. */
    /* 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 (Matrix->DoRealDirect == NULL) {
	if (( Matrix->DoRealDirect = ALLOC(int, Size+1)) == NULL)
	    Matrix->Error = spNO_MEMORY;
    }
    if (Matrix->DoCmplxDirect == NULL) {
	if (( Matrix->DoCmplxDirect = ALLOC(int, Size+1)) == NULL)
	    Matrix->Error = spNO_MEMORY;
    }

    /* Create Intermediate vectors for use in MatrixSolve. */
    if (Matrix->Intermediate == NULL) {
	if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL)