// -*- c++ -*-
// LAPACK++ 1.1 Linear Algebra Package 1.1
// University of Tennessee, Knoxvilee, TN.
// Oak Ridge National Laboratory, Oak Ridge, TN.
// Authors: J. J. Dongarra, E. Greaser, R. Pozo, D. Walker
// (C) 1992-1996 All Rights Reserved
//
// NOTICE
//
// Permission to use, copy, modify, and distribute this software and
// its documentation for any purpose and without fee is hereby granted
// provided that the above copyright notice appear in all copies and
// that both the copyright notice and this permission notice appear in
// supporting documentation.
//
// Neither the Institutions (University of Tennessee, and Oak Ridge National
// Laboratory) nor the Authors make any representations about the suitability
// of this software for any purpose. This software is provided ``as is''
// without express or implied warranty.
//
// LAPACK++ was funded in part by the U.S. Department of Energy, the
// National Science Foundation and the State of Tennessee.
//
// Modifications Copyright (C) 2000-2000, 2002 the R Development Core Team
#ifndef _LA_GEN_FACT_DOUBLE_H
#define _LA_GEN_FACT_DOUBLE_H
#include "lafnames.h"
#include LA_VECTOR_INT_H
#include LA_UNIT_LOWER_TRIANG_MAT_DOUBLE_H
#include LA_UPPER_TRIANG_MAT_DOUBLE_H
#include LA_UTIL_H
#include "lapack.h"
#include "factor.h"
class LaLUFactorDouble : public Factor
{
LaGenMatDouble A_;
LaVectorInt pivot_;
bool singular_;
public:
// constructor
LaLUFactorDouble()
: A_(), pivot_(), singular_(true) { }
inline explicit LaLUFactorDouble(const LaGenMatDouble&);
inline LaLUFactorDouble(const LaLUFactorDouble&);
virtual ~LaLUFactorDouble() { };
// extractor methods for components
LaUnitLowerTriangMatDouble L()
{ LaUnitLowerTriangMatDouble Lmat; Lmat.ref(A_); return Lmat; };
LaUpperTriangMatDouble U()
{ LaUpperTriangMatDouble Umat; Umat.ref(A_); return Umat; };
LaVectorInt pivot()
{ return pivot_; };
bool isSingular()
{ return singular_; };
// linear equation solvers
inline LaGenMatDouble* solve() const;// inverse
inline LaMatDouble& solve(LaMatDouble& B) const; // in-place solution
inline LaMatDouble& solve(LaMatDouble& X, const LaMatDouble& B) const;
// operators
inline LaLUFactorDouble& ref(const LaLUFactorDouble&);
inline LaLUFactorDouble& ref(const LaGenMatDouble&);
};
// constructor/destructor functions
inline LaLUFactorDouble::LaLUFactorDouble(const LaGenMatDouble& A)
: A_(), pivot_(min(A.size(0), A.size(1))), singular_(true)
{
LaGenMatDouble A1;
A1.copy(A);
ref(A1);
}
inline LaLUFactorDouble::LaLUFactorDouble(const LaLUFactorDouble& F)
: A_(), pivot_(), singular_(true)
{
A_.ref(F.A_);
pivot_.ref(F.pivot_);
singular_ = F.singular_;
}
// operators
inline LaLUFactorDouble& LaLUFactorDouble::ref(const LaGenMatDouble& A)
{
if (A.size(0) != A.size(1))
throw(LaException("LaLUFactorDouble::ref(const LaGenMatDouble&)",
"non-square input matrix"));
if(A.inc(0) != 1 || A.inc(1) != 1)
throw(LaException("LaLUFactorDouble::ref(const LaGenMatDouble&)",
"input matrix has non unit increment"));
pivot_.resize(A.size(0));
A_.ref(A);
int info;
F77_CALL(dgetrf)(A.size(0), A.size(0), &A_(0, 0), A.gdim(0),
&pivot_(0), info);
if (info < 0)
throw(LaException("LaLUFactorDouble::ref(const LaGenMatDouble&)",
"illegal input"));
singular_ = info > 0;
return *this;
}
inline LaLUFactorDouble& LaLUFactorDouble::ref(const LaLUFactorDouble& F)
{
A_.ref(F.A_);
pivot_.ref(F.pivot_);
singular_ = F.singular_;
return *this;
}
inline LaGenMatDouble* LaLUFactorDouble::solve() const
{
if (singular_)
throw(LaException("singular matrix"));
int info;
LaGenMatDouble* ans = new LaGenMatDouble();
ans->copy(A_);
int lwork = A_.size(0) *
F77_NAME(ilaenv)(1, "DGETRI", "", A_.size(0), -1, -1, -1);
VectorDouble work(lwork);
F77_CALL(dgetri)(A_.size(0), ans->addr(), A_.gdim(0), pivot_.addr(),
work.addr(), lwork, info);
if (info < 0)
throw(LaException("LaLUFactorDouble::solve()",
"illegal input"));
if (info > 0)
throw(LaException("LaLUFactorDouble::solve()",
"singular matrix"));
return ans;
}
inline LaMatDouble& LaLUFactorDouble::solve(LaMatDouble& B) const
{
if (singular_)
throw(LaException("singular matrix"));
dynamic_cast<LaGenMatDouble&>(B);
int info;
F77_CALL(dgetrs)('N', A_.size(0), B.size(1), A_.addr(),
A_.gdim(0), pivot_.addr(), B.addr(), B.size(0), info);
return B;
}
inline LaMatDouble& LaLUFactorDouble::solve(LaMatDouble& X, const LaMatDouble& B ) const
{
// dynamic_cast<LaGenMatDouble&>(B);
dynamic_cast<LaGenMatDouble&>(X);
X.inject(B);
return solve(X);
}
#if 0
inline void LaGenMatFactorize(LaGenMatDouble &GM, LaLUFactorDouble &GF)
{
integer m = GM.size(0), n = GM.size(1), lda = GM.gdim(0);
integer info=0;
F77NAME(dgetrf)(&m, &n, GM.addr(), &lda, GF.pivot().addr(), &info);
}
inline void LaGenMatFactorizeUnblocked(LaGenMatDouble &A, LaLUFactorDouble &F)
{
integer m = A.size(0), n=A.size(1), lda = A.gdim(0);
integer info=0;
F77NAME(dgetf2)(&m, &n, A.addr(), &lda, F.pivot().addr(), &info);
}
#endif
#endif