C++ GeneralMatrix::numCols方法代码示例

void SymSchurDecomp::getFactor(GeneralMatrix& f) const
{
	if (f.numRows() != q.numRows())
		throw SYLV_MES_EXCEPTION("Wrong dimension of factor matrix in SymSchurDecomp::getFactor");
	if (f.numRows() != f.numCols())
		throw SYLV_MES_EXCEPTION("Factor matrix is not square in SymSchurDecomp::getFactor");
	if (! isPositiveSemidefinite())
		throw SYLV_MES_EXCEPTION("Symmetric decomposition not positive semidefinite in SymSchurDecomp::getFactor");

	f = q;
	for (int i = 0; i < f.numCols(); i++) {
		Vector fi(f, i);
		fi.mult(std::sqrt(lambda[i]));
	}
}

SymSchurDecomp::SymSchurDecomp(const GeneralMatrix& mata)
	: lambda(mata.numRows()), q(mata.numRows())
{
	// check mata is square
	if (mata.numRows() != mata.numCols())
		throw SYLV_MES_EXCEPTION("Matrix is not square in SymSchurDecomp constructor");

	// prepare for dsyevr
	const char* jobz = "V";
	const char* range = "A";
	const char* uplo = "U";
	lapack_int n = mata.numRows();
	GeneralMatrix tmpa(mata);
	double* a = tmpa.base();
	lapack_int lda = tmpa.getLD();
	double dum;
	double* vl = &dum;
	double* vu = &dum;
	lapack_int idum;
	lapack_int* il = &idum;
	lapack_int* iu = &idum;
	double abstol = 0.0;
	lapack_int m = n;
	double* w = lambda.base();
	double* z = q.base();
	lapack_int ldz = q.getLD();
	lapack_int* isuppz = new lapack_int[2*std::max(1,(int) m)];
	double tmpwork;
	lapack_int lwork = -1;
	lapack_int tmpiwork;
	lapack_int liwork = -1;
	lapack_int info;

	// query for lwork and liwork
	dsyevr(jobz, range, uplo, &n, a, &lda, vl, vu, il, iu, &abstol,
				  &m, w, z, &ldz, isuppz, &tmpwork, &lwork, &tmpiwork, &liwork, &info);
	lwork = (int)tmpwork;
	liwork = tmpiwork;
	// allocate work arrays
	double* work = new double[lwork];
	lapack_int* iwork = new lapack_int[liwork];
	
	// do the calculation
	dsyevr(jobz, range, uplo, &n, a, &lda, vl, vu, il, iu, &abstol,
				  &m, w, z, &ldz, isuppz, work, &lwork, iwork, &liwork, &info);

	if (info < 0)
		throw SYLV_MES_EXCEPTION("Internal error in SymSchurDecomp constructor");
	if (info > 0)
		throw SYLV_MES_EXCEPTION("Internal LAPACK error in DSYEVR");

	delete [] work;
	delete [] iwork;
	delete [] isuppz;
}

void SylvMatrix::multLeft(int zero_cols, const GeneralMatrix& a, const GeneralMatrix& b)
{
	int off = a.numRows() - a.numCols();
	if (off < 0 || a.numRows() != rows || off != zero_cols ||
		rows != b.numRows() || cols != b.numCols()) {
		throw SYLV_MES_EXCEPTION("Wrong matrix dimensions for multLeft.");
	}
	// here we cannot call SylvMatrix::gemm since it would require
	// another copy of (usually big) b (we are not able to do inplace
	// submatrix of const GeneralMatrix)
	if (a.getLD() > 0 && ld > 0) {
		blas_int mm = a.numRows();
		blas_int nn = cols;
		blas_int kk = a.numCols();
		double alpha = 1.0;
		blas_int lda = a.getLD();
		blas_int ldb = ld;
		double beta = 0.0;
		blas_int ldc = ld;
		dgemm("N", "N", &mm, &nn, &kk, &alpha, a.getData().base(), &lda,
				   b.getData().base()+off, &ldb, &beta, data.base(), &ldc);
	}
}

void SqSylvMatrix::multInvLeft2(GeneralMatrix& a, GeneralMatrix& b,
								double& rcond1, double& rcondinf) const
{
	if (rows != a.numRows() || rows != b.numRows()) {
		throw SYLV_MES_EXCEPTION("Wrong dimensions for multInvLeft2.");
	}
	// PLU factorization
	Vector inv(data);
	lapack_int * const ipiv = new lapack_int[rows];
	lapack_int info;
	lapack_int rows2 = rows;
	dgetrf(&rows2, &rows2, inv.base(), &rows2, ipiv, &info);
	// solve a
	lapack_int acols = a.numCols();
	double* abase = a.base();
	dgetrs("N", &rows2, &acols, inv.base(), &rows2, ipiv,
				  abase, &rows2, &info);
	// solve b
	lapack_int bcols = b.numCols();
	double* bbase = b.base();
	dgetrs("N", &rows2, &bcols, inv.base(), &rows2, ipiv,
				  bbase, &rows2, &info);
	delete [] ipiv;

	// condition numbers
	double* const work = new double[4*rows];
	lapack_int* const iwork = new lapack_int[rows];
	double norm1 = getNorm1();
	dgecon("1", &rows2, inv.base(), &rows2, &norm1, &rcond1, 
				  work, iwork, &info);
	double norminf = getNormInf();
	dgecon("I", &rows2, inv.base(), &rows2, &norminf, &rcondinf, 
				  work, iwork, &info);
	delete [] iwork;
	delete [] work;
}

Vector::Vector(int row, GeneralMatrix& m)
	: len(m.numCols()), s(m.getLD()), data(&(m.get(row, 0))), destroy(false)
{
}

SchurDecompZero::SchurDecompZero(const GeneralMatrix& m)
	: SchurDecomp(SqSylvMatrix(m, m.numRows()-m.numCols(), 0, m.numCols())),
	  ru(m, 0, 0, m.numRows()-m.numCols(), m.numCols())
{
	ru.multRight(getQ());
}