C++ DMat::data方法代码示例

//---------------------------------------------------------
void eig(const DMat& A, DVec& Re, DMat& VL, DMat& VR, bool bL, bool bR)
//---------------------------------------------------------
{
  // Compute eigensystem of a real general matrix
  // Currently NOT returning imaginary components

  static DMat B;

  if (!A.is_square()) { umERROR("eig(A)", "matrix is not square."); }

  int N = A.num_rows();
  int LDA=N, LDVL=N, LDVR=N, ldwork=10*N, info=0;

  Re.resize(N);     // store REAL components of eigenvalues in Re
  VL.resize(N,N);   // storage for LEFT eigenvectors
  VR.resize(N,N);   // storage for RIGHT eigenvectors
  DVec Im(N);     // NOT returning imaginary components
  DVec work(ldwork, 0.0);

  // Work on a copy of A
  B = A;

  char jobL = bL ? 'V' : 'N';   // calc LEFT eigenvectors?
  char jobR = bR ? 'V' : 'N';   // calc RIGHT eigenvectors?

  GEEV (jobL,jobR, N, B.data(), LDA, Re.data(), Im.data(), 
        VL.data(), LDVL, VR.data(), LDVR, work.data(), ldwork, info);

  if (info < 0) { 
    umERROR("eig(A, Re,Im)", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    umLOG(1, "eig(A, Re,Im): ...\n"
             "\nThe QR algorithm failed to compute all the"
             "\neigenvalues, and no eigenvectors have been" 
             "\ncomputed;  elements %d+1:N of WR and WI contain"
             "\neigenvalues which have converged.\n", info);
  }

#if (0)
  // Return (Re,Imag) parts of eigenvalues as columns of Ev
  Ev.resize(N,2);
  Ev.set_col(1, Re);
  Ev.set_col(2, Im);
#endif

#ifdef _DEBUG
    //#####################################################
    // Check for imaginary components in eigenvalues
    //#####################################################
    double im_max = Im.max_val_abs();
    if (im_max > 1e-6) {
      umERROR("eig(A)", "imaginary components in eigenvalues.");
    }
    //#####################################################
#endif
}

//---------------------------------------------------------
DMat& qr(DMat& A, bool in_place)
//---------------------------------------------------------
{
  // Form orthogonal QR factorization of A(m,n). 
  // The result Q is represented as a product of 
  // min(m, n) elementary reflectors. 

  int M=A.num_rows(), N=A.num_cols(), LDA=A.num_rows();
  int min_mn = A.min_mn(), info=0; DVec tau(min_mn);

  if (in_place) 
  {
    // factorize arg
    GEQRF(M, N, A.data(), LDA, tau.data(), info);

    if (info) { umERROR("qr(A)", "dgeqrf reports: info = %d", info); }
  //A.set_qrtau(tau);         // H(i) = I - tau * v * v'
    A.set_factmode(FACT_QR);  // indicate factored state
    return A;
  } 
  else
  {
    // factorize copy of arg
    DMat* tmp = new DMat(A, OBJ_temp, "qr(A)");
    GEQRF (M, N, tmp->data(), LDA, tau.data(), info);

    if (info) { umERROR("qr(A)", "dgeqrf reports: info = %d", info); }
  //tmp->set_qrtau(tau);         // H(i) = I - tau * v * v'
    tmp->set_factmode(FACT_QR);  // indicate factored state
    return (*tmp);
  }
}

//---------------------------------------------------------
DMat& chol(DMat& A, bool in_place)
//---------------------------------------------------------
{
  // Given symmetric positive-definite matrix A,
  // return its Cholesky-factorization for use
  // later in solving (multiple) linear systems.

  int M=A.num_rows(), LDA=A.num_rows(), info=0;
  char uplo = 'U';

  if (in_place) 
  {
    // factorize arg
    POTRF (uplo, M, A.data(), LDA, info);
    if (info) { umERROR("chol(A)", "dpotrf reports: info = %d", info); }
    A.zero_below_diag();
    A.set_factmode(FACT_CHOL);  // indicate factored state
    return A;
  } 
  else
  {
    // factorize copy of arg
    DMat* tmp = new DMat(A, OBJ_temp, "chol(A)");
    POTRF (uplo, M, tmp->data(), LDA, info);
    if (info) { umERROR("chol(A)", "dpotrf reports: info = %d", info); }
    tmp->zero_below_diag();
    tmp->set_factmode(FACT_CHOL);  // indicate factored state
#if (0)
    // compare with Matlab
    tmp->print(g_MSGFile, "chol", "lf", 4, 8);
#endif
    return (*tmp);
  }
}

//---------------------------------------------------------
DMat& lu(DMat& A, bool in_place)
//---------------------------------------------------------
{
  // Given square matrix A, return its lu-factorization 
  // for use later in solving (multiple) linear systems.

  if (!A.is_square()) { umERROR("lu(A)", "matrix is not square."); }
  int rows=A.num_rows(); int N=rows, LDA=rows, info=0;
  int* ipiv = umIVector(rows);

  if (in_place) 
  {
    // factorize arg
    GETRF(N, N, A.data(), LDA, ipiv, info);
    if (info) { umERROR("lu(A)", "dgetrf reports: info = %d", info); }
    A.set_pivots(ipiv);        // store pivots
    A.set_factmode(FACT_LUP);  // indicate factored state
    return A;
  } 
  else
  {
    // factorize copy of arg
    DMat* tmp = new DMat(A, OBJ_temp, "lu(A)");
    GETRF(N, N, tmp->data(), LDA, ipiv, info);
    if (info) { umERROR("lu(A)", "dgetrf reports: info = %d", info); }
    tmp->set_pivots(ipiv);        // store pivots
    tmp->set_factmode(FACT_LUP);  // indicate factored state
    return (*tmp);
  }
}

// compute eigensystem of a real symmetric matrix
//---------------------------------------------------------
void eig_sym(const DMat& A, DVec& ev, DMat& Q, bool bDoEVecs)
//---------------------------------------------------------
{
  if (!A.is_square()) { umERROR("eig_sym(A)", "matrix is not square."); }

  int N = A.num_rows();
  int LDA=N, LDVL=N, LDVR=N, ldwork=10*N, info=0;
  DVec work(ldwork, 0.0, OBJ_temp, "work_TMP");

  Q = A;          // Calculate eigenvectors in Q (optional)
  ev.resize(N);   // Calculate eigenvalues in ev

  char jobV = bDoEVecs ? 'V' : 'N';

  SYEV (jobV,'U', N, Q.data(), LDA, ev.data(), work.data(), ldwork, info);  

  if (info < 0) { 
    umERROR("eig_sym(A, Re,Im)", "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    umLOG(1, "eig_sym(A, W): ...\n"
             "\nthe algorithm failed to converge;"
             "\n%d off-diagonal elements of an intermediate"
             "\ntridiagonal form did not converge to zero.\n", info);
  }
}

// DPOSV uses Cholesky factorization A=U^T*U, A=L*L^T 
// to compute the solution to a real system of linear 
// equations A*X=B, where A is a square, (N,N) symmetric 
// positive definite matrix and X and B are (N,NRHS).
//
// If the system is over or under-determined, 
// (i.e. A is not square), then pass the problem
// to the Least-squares solver (DGELSS) below.
//---------------------------------------------------------
void umSOLVE_CH(const DMat& mat, const DMat& B, DMat& X)
//---------------------------------------------------------
{
  if (!mat.ok()) {umWARNING("umSOLVE_CH()", "system is empty"); return;}
  if (!mat.is_square()) {
    umSOLVE_LS(mat, B, X);    // return a least-squares solution.
    return;
  }
  
  DMat A(mat);    // Work with a copy of input array.
  X = B;          // initialize solution with rhs

  int rows=A.num_rows(), LDA=A.num_rows(), cols=A.num_cols();
  int LDB=X.num_rows(), NRHS=X.num_cols(), info=0;
  assert(LDB >= rows);  // enough space for solutions?

  // Solve the system.
  POSV('U', rows, NRHS, A.data(), LDA, X.data(), LDB, info);

  if (info < 0) { 
    X = 0.0;
    umERROR("umSOLVE_CH(A,B, X)", 
            "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    X = 0.0;
    umERROR("umSOLVE_CH(A,B, X)", 
            "\nINFO = %d.  The leading minor of order %d of A"
            "\nis not positive definite, so the factorization" 
            "\ncould not be completed. No solution computed.", 
              info, info);
  }
}

// DGESV uses the LU factorization to compute solution 
// to a real system of linear equations, A * X = B, 
// where A is square (N,N) and X, B are (N,NRHS).
//
// If the system is over or under-determined, 
// (i.e. A is not square), then pass the problem
// to the Least-squares solver (DGELSS) below.
//---------------------------------------------------------
void umSOLVE(const DMat& mat, const DMat& B, DMat& X)
//---------------------------------------------------------
{
  if (!mat.ok()) {umWARNING("umSOLVE()", "system is empty"); return;}
  if (!mat.is_square()) {
    umSOLVE_LS(mat, B, X);    // return a least-squares solution.
    return;
  }

  DMat A(mat);    // work with copy of input
  X = B;          // initialize result with RHS

  int rows=A.num_rows(), LDA=A.num_rows(), cols=A.num_cols();
  int LDB=B.num_rows(), NRHS=B.num_cols(), info=0;
  if (rows<1) {umWARNING("umSOLVE()", "system is empty"); return;}
  IVec ipiv(rows);

  // Solve the system.
  GESV(rows, NRHS, A.data(), LDA, ipiv.data(), X.data(), LDB, info);

  if (info < 0) { 
    X = 0.0;
    umERROR("umSOLVE(A,B, X)", 
            "Error in input argument (%d)\nNo solution computed.", -info);
  } else if (info > 0) {
    X = 0.0;
    umERROR("umSOLVE(A,B, X)", 
            "\nINFO = %d.  U(%d,%d) was exactly zero."
            "\nThe factorization has been completed, but the factor U is "
            "\nexactly singular, so the solution could not be computed.", 
              info, info, info);
  }
}

//---------------------------------------------------------
bool chol_solve(const DMat& ch, const DMat& B, DMat& X)
//---------------------------------------------------------
{
  // Solve a set of linear systems using Cholesky-factored 
  // symmetric positive-definite matrix, A = U^T U.

  if (FACT_CHOL != ch.get_factmode()) {umERROR("chol_solve(ch,B,X)", "matrix is not factored.");}
  int M =ch.num_rows(), lda=ch.num_rows(); 
  int ldb=B.num_rows(), nrhs=B.num_cols(); assert(ldb == M);
  char uplo = 'U';  int info=0; 
  double* ch_data = const_cast<double*>(ch.data());

  X = B;  // overwrite X with RHS's, then solutions
  POTRS (uplo, M, nrhs, ch_data, lda, X.data(), ldb, info);

  if (info) { umERROR("chol_solve(ch,B,X)", "dpotrs reports: info = %d", info); }
  return true;
}

//---------------------------------------------------------
void umAxB(const DMat& A, const DMat& B, DMat& C)
//---------------------------------------------------------
{
  //-------------------------
  // C = A * B
  //-------------------------
  // A = op(A) is (M,K)
  // B = op(B) is (K,N)
  //        C  is (M,N)
  //-------------------------
  int M=A.num_rows(), K=A.num_cols(), N=B.num_cols();
  int LDA=M, LDB=K, LDC=M;
  double one=1.0, zero=0.0;
  if (B.num_rows() != K) { umERROR("umAxB(A,B,C)", "wrong dimensions"); }
  C.resize(M,N);

  GEMM ('N','N',M,N,K, one,A.data(),LDA, 
                           B.data(),LDB, 
                      zero,C.data(),LDC);
}

// DGELSS computes minimum norm solution to a real linear 
// least squares problem:   Minimize 2-norm(| b - A*x |).   
// using the singular value decomposition (SVD) of A. 
// A is an M-by-N matrix which may be rank-deficient.   
//---------------------------------------------------------
void umSOLVE_LS(const DMat& mat, const DMat& B, DMat& X)
//---------------------------------------------------------
{
  if (!mat.ok()) {umWARNING("umSOLVE_LS()", "system is empty"); return;}

  DMat A(mat);    // work with copy of input.

  int rows=A.num_rows(), cols=A.num_cols(), mmn=A.min_mn();
  int LDB=A.max_mn(), NRHS=B.num_cols();
  if (rows!=B.num_rows()) {umERROR("umSOLVE_LS(A,B)", "Inconsistant matrix sizes.");}

  DVec s(mmn);    // allocate array for singular values

  // X must be big enough to store various results.
  // Resize X so that its leading dimension = max(M,N), 
  // then load the set of right hand sides.

  X.resize(LDB,NRHS, true, 0.0);

  for (int j=1; j<=NRHS; ++j)     // loop across colums
    for (int i=1; i<=rows; ++i)   // loop down rows
      X(i,j) = B(i,j);

  // RCOND is used to determine the effective rank of A.   
  // Singular values S(i) <= RCOND*S(1) are treated as zero.   
  // If RCOND < 0, machine precision is used instead.   

//double rcond =  1.0 / 1.0e16;
  double rcond = -1.0;

  // NBN: ACML does not use the work vector.
  int mnLo=A.min_mn(), mnHi=A.max_mn(), rank=1, info=1;
  int lwork = 10*mnLo + std::max(2*mnLo, std::max(mnHi, NRHS));
  DVec work(lwork); 

  // Solve the system
  GELSS (rows, cols, NRHS, A.data(), rows, X.data(), LDB, s.data(), rcond, rank, work.data(), lwork, info);

  //---------------------------------------------
  // Report:
  //---------------------------------------------

  if (info == 0) {
    umLOG(1, "umSOLVE_LS reports successful LS-solution."
             "\nRCOND = %0.6e, "
             "\nOptimal length of work array was %d\n", rcond, lwork);
  } 
  else 
  {
    if (info < 0) { 
      X = 0.0;
      umERROR("umSOLVE_LS(DMat&, DMat&)", 
              "Error in input argument (%d)\nNo solution or error bounds computed.", -info);

    } else if (info > 0) {
      X = 0.0;
      umERROR("umSOLVE_LS(DMat&, DMat&)", 
          "\nThe algorithm for computing the SVD failed to converge.\n"
          "\n%d off-diagonal elements of an intermediate "
          "\nbidiagonal form did not converge to zero.\n "
          "\nRCOND = %0.6e, "
          "\nOptimal length of work array was %d.\n", info, rcond, lwork);
    }
  }
}