linslv.cc

LAPACK++ (Linear Algebra PACKage in C++) is a software library for numerical linear algebra that sol

}



    

void LaQRLinearSolve(const LaGenMatComplex& A, LaGenMatComplex& X, 
		     const LaGenMatComplex& B )
{
    LaGenMatComplex A1(A);
    LaQRLinearSolveIP(A1, X, B);
}

    
// General QR solver
//
//                          M x N              N x nrhs           M  x nrhs
//
void LaQRLinearSolveIP(LaGenMatComplex& A, LaGenMatComplex& X, const LaGenMatComplex& B )
{
#ifndef HPPA
    const char fname[] = "LaQRLinearSolveIP(LaGenMatComplex &A, &X, &B)";
#else
    char *fname = NULL;  // HP C++ does not support string initalization!
#endif

    // let's not worry about non-unit column strides for the moment
    if ( A.inc(0) != 1 || A.inc(1) != 1)
        throw(LaException(fname, "A is  non-contiguous."));
    if ( A.size(0) == 0 || A.size(1) == 0 )
        throw(LaException(fname, "Matrix A is empty; one dimension is zero."));

    if (!(  A.size(0) == B.size(0) &&
            A.size(1) == X.size(0) &&
            X.size(1) == B.size(1) ))
        throw(LaException(fname, "input matrices are non-conformant."));

    long int info = 0;
    int M = A.size(0);
    int N = A.size(1);
    long int Ml = M;
    long int Nl = N;
    int nrhs = X.size(1);
    long int nrhsl = nrhs;
    long int lda = A.inc(0) * A.gdim(0);

        
    //int nb = 32;
    int nb = LaEnvBlockSize("ZGELS", A);
    long int lwork = M * N + nb * std::max(M*N, nrhs);
    //std::cout << "Block size: " << nb << std::endl;

    LaVectorComplex WORK(lwork);         
    char trans = 'N';
    
    if (M != N)
    {
	// Typically is A non-square, so we need to create tmp X because 
	// X is N x nrhs, while B is M x nrhs.  We need to make copies of
	// these so that the routine won't corrupt data around X and B.

        LaGenMatComplex Xtmp(std::max(M, N), nrhs);
        long int ldx = Xtmp.inc(0) * Xtmp.gdim(0);

	// Copy B into the temporary X matrix that is passed to zgels()
        Xtmp(LaIndex(0,M-1),LaIndex()).inject( B );

        F77NAME(zgels) (&trans, &Ml, &Nl, &nrhsl, &A(0,0), &lda, &Xtmp(0,0), 
                &ldx, &WORK(0), &lwork, &info);

	// And copy the result from the larger matrix back into
	// the actual result matrix.
	X.inject(Xtmp(LaIndex(0,N-1), LaIndex()));
    }
    else
    {
        long int ldx = X.inc(0) * X.gdim(0);

	// Copy B into the X matrix that is passed to dgels()
        X.inject( B );

        F77NAME(zgels) (&trans, &Ml, &Nl, &nrhsl, &A(0,0), &lda, &X(0,0), 
                &ldx, &WORK(0), &lwork, &info);
    }


    // this shouldn't really happen.
    //
    if (info < 0)
        throw(LaException(fname, "Internal error in LAPACK: ZGELS()"));

}

#endif // LA_COMPLEX_SUPPORT



// ////////////////////////////////////////////////////////////

void  LaCholLinearSolve( const LaSpdMatDouble& A, LaGenMatDouble& X,
        LaGenMatDouble& B )
{
        LaSpdMatDouble A1(A);
        LaCholLinearSolveIP(A1, X, B);
}

void  LaCholLinearSolve( const LaSymmMatDouble& A, LaGenMatDouble& X,
			 const LaGenMatDouble& B )
{
        LaSymmMatDouble A1(A);
        LaCholLinearSolveIP(A1, X, B);
}

// A is NxN, X is MxN and B is MxN
//
void LaCholLinearSolveIP( LaSpdMatDouble& A, LaGenMatDouble& X, 
        LaGenMatDouble& B )
{
#ifndef HPPA
    const char fname[] = "LaCholLinearSolveIP(LaSpdMatDouble &A, &X, &B)";
#else
    char *fname = NULL;  // HP C++ does not support string initalization!
#endif

    // let's not worry about non-unit column strides for the moment
    if ( A.inc(0) != 1 || A.inc(1) != 1)
        throw(LaException(fname, "A is  non-contiguous."));

    if (!(X.size(0) == B.size(0) && X.size(1) == B.size(1)))
        throw(LaException(fname, "X and B are non-conformant."));

    X.inject(B);            // will throw exception if not conformant


    // in the future this can call the linear least square routines
    // to handle non-square matrices

    if (A.size(0) != A.size(1))
        throw(LaException(fname, "Square matrix expected.\n"));

    if (A.size(1) != X.size(0))
        throw(LaException(fname, "A and X are non-comformant."));

    long int info = 0;
    long int M = A.size(0);
    //long int N = A.size(1);
    long int K = X.size(1);
    long int lda = A.inc(0) * A.gdim(0);
    long int ldx = X.inc(0) * X.gdim(0);
    char uplo = 'L';

    F77NAME(dposv) (&uplo, &M, &K, &A(0,0), &lda,  &X(0,0), &ldx, &info);

    // this shouldn't really happen.
    //
    if (info < 0)
        throw(LaException(fname, "Internal error in LAPACK: SGESV()"));

    if (info > 0)
        throw (LaException(fname, "A is not symmetric-positive-definite."));


}

void LaCholLinearSolveIP( LaSymmMatDouble& A, LaGenMatDouble& X, 
			  const LaGenMatDouble& B )
{
#ifndef HPPA
    const char fname[] = "LaCholLinearSolveIP(LaSymmMatDouble &A, &X, &B)";
#else
    char *fname = NULL;  // HP C++ does not support string initalization!
#endif

    // let's not worry about non-unit column strides for the moment
    if ( A.inc(0) != 1 || A.inc(1) != 1)
        throw(LaException(fname, "A is  non-contiguous."));

    if (!(X.size(0) == B.size(0) && X.size(1) == B.size(1)))
        throw(LaException(fname, "X and B are non-conformant."));

    X.inject(B);            // will throw exception if not conformant


    // in the future this can call the linear least square routines
    // to handle non-square matrices

    if (A.size(0) != A.size(1))
        throw(LaException(fname, "Square matrix expected.\n"));

    if (A.size(1) != X.size(0))
        throw(LaException(fname, "A and X are non-comformant."));

    long int info = 0;
    integer M = A.size(0);
    //integer N = A.size(1);
    integer K = X.size(1);
    integer lda = A.inc(0) * A.gdim(0);
    integer ldx = X.inc(0) * X.gdim(0);
    char uplo = 'L';

    LaVectorLongInt ipiv(M);
    integer lwork = -1;
    LaVectorDouble work(1);
    // Workspace query
    F77NAME(dsysv) (&uplo, &M, &K, &A(0,0), &lda,  &ipiv(0), &X(0,0), &ldx, 
		    &work(0), &lwork, &info);
    lwork = integer(work(0));
    work.resize(lwork, 1);

    F77NAME(dsysv) (&uplo, &M, &K, &A(0,0), &lda,  &ipiv(0), &X(0,0), &ldx, 
		    &work(0), &lwork, &info);

    // this shouldn't really happen.
    //
    if (info < 0)
        throw(LaException(fname, "Internal error in LAPACK: DSYSV()"));

    if (info > 0)
        throw (LaException(fname, "Matrix is singular."));


}


void LUFactorizeIP(LaGenMatDouble &GM, LaVectorLongInt &PIV)
{
    integer m = GM.size(0), n = GM.size(1), lda = GM.gdim(0);
    integer info=0;
    assert(PIV.size() >= (m<n ? m : n));

    // Copied from LaGenMatFactorize in fmd.h of LAPACK++, but the
    // version there didn't admit to being in-place modifying.
    F77NAME(dgetrf)(&m, &n, &GM(0,0), &lda, &(PIV(0)), &info);

    if (info < 0)
	throw LaException("LUFactorizeIP", "Error in argument");
}

#ifdef LA_COMPLEX_SUPPORT
void LUFactorizeIP(LaGenMatComplex &GM, LaVectorLongInt &PIV)
{
    integer m = GM.size(0), n = GM.size(1), lda = GM.gdim(0);
    integer info=0;
    assert(PIV.size() >= (m<n ? m : n));

    // Copied from LaGenMatFactorize in fmd.h of LAPACK++, but the
    // version there didn't admit to being in-place modifying.
    F77NAME(zgetrf)(&m, &n, &GM(0,0), &lda, &(PIV(0)), &info);

    if (info < 0)
	throw LaException("LUFactorizeIP", "Error in argument");
}

void LaLUInverseIP(LaGenMatComplex &A, LaVectorLongInt &PIV)
{
    LaVectorComplex work; // will be resized in other function
    LaLUInverseIP(A, PIV, work);
}
void LaLUInverseIP(LaGenMatComplex &A, LaVectorLongInt &PIV, LaVectorComplex &work)
{
    integer N = A.size(1), lda = A.gdim(0), info = 0;
    if(A.size(0) != A.size(1))
      throw LaException("LaLUInverseIP", "Input must be square");
    long int W = work.size();

    // Check for minimum work size
    if ( W < A.size(0) ) 
    {
      int nb = LaEnvBlockSize("ZGETRI", A);
      W = N*nb;
      work.resize(W, 1);
    }

    F77NAME(zgetri)(&N, &(A(0,0)), &lda, &(PIV(0)), &work(0), &W, &info);
    if (info < 0)
	throw LaException("LaLUInverseIP", "Error in zgetri argument");
    if (info > 0)
	throw LaException("LaLuInverseIP", "Matrix is singlular, cannot compute inverse");
}
#endif // LA_COMPLEX_SUPPORT

void LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV)
{
    LaVectorDouble work; // will be resized in other function
    LaLUInverseIP(A, PIV, work);
#if 0
    // above code is shorter - remove this code soon.
    integer N = A.size(1), lda = A.gdim(0), info = 0;
    if(A.size(0) != A.size(1))
      throw LaException("LaLUInverseIP", "Input must be square");
    int nb = LaEnvBlockSize("DGETRI", A);
    long int W = N*nb;
    LaVectorDouble work(W);         

    F77NAME(dgetri)(&N, &(A(0,0)), &lda, &(PIV(0)), &work(0), &W, &info);
    if (info < 0)
	throw LaException("inv", "Error in dgetri argument");
#endif
}

void LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, LaVectorDouble &work)
{
    integer N = A.size(1), lda = A.gdim(0), info = 0;
    if(A.size(0) != A.size(1))
      throw LaException("LaLUInverseIP", "Input must be square");
    long int W = work.size();

    // Check for minimum work size
    if ( W < A.size(0) ) 
    {
      int nb = LaEnvBlockSize("DGETRI", A);
      W = N*nb;
      work.resize(W, 1);
    }

    F77NAME(dgetri)(&N, &(A(0,0)), &lda, &(PIV(0)), &work(0), &W, &info);
    if (info < 0)
	throw LaException("LaLUInverseIP", "Error in dgetri argument");
    if (info > 0)
	throw LaException("LaLuInverseIP", "Matrix is singlular, cannot compute inverse");
}