spfactor.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

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






/*!
 *  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 \a Mode.  If \a Mode is set
 *  to \a spDIRECT_PARTITION, then the all rows are placed in the direct
 *  addressing partition.  Similarly, if \a Mode is set to
 *  \a spINDIRECT_PARTITION, then the all rows are placed in the indirect
 *  addressing partition.  By setting \a Mode to \a 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.
 *
 *  \param eMatrix
 *      Pointer to matrix.
 *  \param Mode
 *      Mode must be one of three special codes: \a spDIRECT_PARTITION,
 *      \a spINDIRECT_PARTITION, or \a spAUTO_PARTITION.
 */

void
spPartition(
    spMatrix eMatrix,
    int Mode
)
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement, pColumn;
register  int  Step, Size;
register  int  *Nc, *No;
register  long  *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 vASSERT( Mode == spAUTO_PARTITION, "Invalid partition code" );;

/* Otherwise, count all operations needed in when factoring matrix. */
    Nc = (int *)Matrix->MarkowitzRow;
    No = (int *)Matrix->MarkowitzCol;
    Nm = (long *)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.
 *
 *  >>> Possible errors: