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

// 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).
//---------------------------------------------------------
void umSOLVE_CH(const DMat& mat, const DVec& b, DVec& x)
//---------------------------------------------------------
{
  // check args
  assert(mat.is_square());            // symmetric
  assert(b.size() >= mat.num_rows()); // is b consistent?
  assert(b.size() <= x.size());       // can x store solution?
  
  DMat A(mat);    // work with copy of input
  x = b;          // allocate solution vector

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

  // 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);
  }
}

//---------------------------------------------------------
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
}