C++ DMat类代码示例

//---------------------------------------------------------
void xyztorst
(
  const DVec& X,  // [in]
  const DVec& Y,  // [in]
  const DVec& Z,  // [in]
        DVec& r,  // [out]
        DVec& s,  // [out]
        DVec& t   // [out]
)
//---------------------------------------------------------
{
  // function [r,s,t] = xyztorst(x, y, z)
  // Purpose : Transfer from (x,y,z) in equilateral tetrahedron
  //           to (r,s,t) coordinates in standard tetrahedron

  double sqrt3=sqrt(3.0), sqrt6=sqrt(6.0); int Nc=X.size();
  DVec v1(3),v2(3),v3(3),v4(3);
  DMat tmat1(3,Nc), A(3,3), rhs; 
  
  v1(1)=(-1.0);  v1(2)=(-1.0/sqrt3);  v1(3)=(-1.0/sqrt6);
  v2(1)=( 1.0);  v2(2)=(-1.0/sqrt3);  v2(3)=(-1.0/sqrt6);
  v3(1)=( 0.0);  v3(2)=( 2.0/sqrt3);  v3(3)=(-1.0/sqrt6);
  v4(1)=( 0.0);  v4(2)=( 0.0      );  v4(3)=( 3.0/sqrt6);

  // back out right tet nodes
  tmat1.set_row(1,X); tmat1.set_row(2,Y); tmat1.set_row(3,Z);
  rhs = tmat1 - 0.5*outer(v2+v3+v4-v1, ones(Nc));
  A.set_col(1,0.5*(v2-v1)); A.set_col(2,0.5*(v3-v1)); A.set_col(3,0.5*(v4-v1));

  DMat RST = A|rhs;

  r=RST.get_row(1); s=RST.get_row(2); t=RST.get_row(3);
}

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

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

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

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

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

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

bool isZero(const DMat<RT>& A)
{
  size_t len = A.numRows() * A.numColumns();
  if (len == 0) return true;
  for (auto t = A.array() + len - 1; t >= A.array(); t--)
    if (!A.ring().is_zero(*t)) return false;
  return true;
}

//---------------------------------------------------------
void Euler3D::CouetteBC3D
(
    const DVec& xi,
    const DVec& yi,
    const DVec& zi,
    const DVec& nxi,
    const DVec& nyi,
    const DVec& nzi,
    const IVec& tmapI,
    const IVec& tmapO,
    const IVec& tmapW,
    const IVec& tmapC,
    double ti,
    DMat& Qio
)
//---------------------------------------------------------
{
    //#######################################################
    // FIXME: what about top and bottom of 3D annulus?
    //#######################################################


    // gmapB = concat(mapI,mapO,mapW, ...);
    // Check for no boundary faces
    if (gmapB.size() < 1) {
        return;
    }

    // function [Q] = CouetteBC3D(xin, yin, nxin, nyin, mapI, mapO, mapW, mapC, Q, time);
    // Purpose: evaluate solution for Couette flow

    // Couette flow (mach .2 at inner cylinder)
    // gamma = 1.4;

    int Nr = Qio.num_rows();
    DVec rho,rhou,rhov,rhow,Ener;

    // wrap current state data (columns of Qio)
    rho.borrow (Nr, Qio.pCol(1));
    rhou.borrow(Nr, Qio.pCol(2));
    rhov.borrow(Nr, Qio.pCol(3));
    rhow.borrow(Nr, Qio.pCol(4));
    Ener.borrow(Nr, Qio.pCol(5));

    // update boundary nodes of Qio with boundary data
    // pre-calculated in function precalc_bdry_data()

    rho (gmapB) = rhoB;
    rhou(gmapB) = rhouB;
    rhov(gmapB) = rhovB;
    rhow(gmapB) = rhowB;
    Ener(gmapB) = EnerB;
}

//---------------------------------------------------------
void CurvedINS2D::Report(bool bForce)
//---------------------------------------------------------
{
  static DMat vort;
  bool normal_report=true;
  if (eVolkerCylinder == sim_type) {normal_report=false;}

#if (1)
  normal_report=true;
#endif

  // print report header on first step
  if (1 == tstep) {
    if (normal_report) {
      umLOG(1, "\n  step    time      min(Ux)    max(Ux)   min(Vort)    max(Vort)\n"
                 "--------------------------------------------------------------------\n");
    } else {
      umLOG(1, "\n  step     time        Cd/ra      Cl/ra       dP   \n"
                 "---------------------------------------------------\n");
    }
  }

  if (normal_report) {
    if (!umMOD(tstep,Nreport)||(1==tstep)||(tstep==Nsteps)||bForce) {

      Curl2D(Ux, Uy, vort);

      umLOG(1, "%7d  %6.3lf  %10.5lf %10.5lf  %10.2lf %10.2lf\n", 
               tstep, time, Ux.min_val(), Ux.max_val(), vort.min_val(), vort.max_val());
    }
  } else {
    // VolkerCylinder: compute coefficients of drag and lift,
    // as well as the pressure drop at the two sample points.
    INSLiftDrag2D(0.05);
  }


  if (!umMOD(tstep,Nrender)||(1==tstep)||(tstep==Nsteps)||bForce) 
  {
    // (this->*ExactSolution)(x, y, time, nu, exUx, exUy, exPR);

    Curl2D(Ux, Uy, vort);

    // load render data
    Q_plot.set_col(1, Ux); // -exUx);
    Q_plot.set_col(2, Uy); // -exUy);
    Q_plot.set_col(3, PR); // -exPR);
    Q_plot.set_col(4, vort);

    OutputVTK(Q_plot, NvtkInterp);
  }
}

void addTo(DMat<RT>& A, MatrixWindow wA, const DMat<RT>& B, MatrixWindow wB)
{
  assert(wA.sameSize(wB));
  long rA = wA.begin_row;
  long rB = wB.begin_row;
  for (; rA < wA.end_row; ++rA, ++rB)
    {
      long cA = wA.begin_column;
      long cB = wB.begin_column;
      for (; cA < wA.end_column; ++cA, ++cB)
        {
          auto& a = A.entry(rA, cA);
          A.ring().add(a, a, B.entry(rB, cB));
        }
    }
}

//---------------------------------------------------------
void CurvedCNS2D::BoxFlowIC2D
(
  const DVec&   xi,   // [in]
  const DVec&   yi,   // [in]
        double  ti,   // [in]
        DMat&   Qo    // [out]
)
//---------------------------------------------------------
{
  // function Q = BoxFlowIC2D(x, y, time)
  //
  // Purpose: compute plane flow configuration 

  //-------------------------------------
  // adjust parameters
  //-------------------------------------
  this->gamma = 1.4;
  this->gm1   = 0.4;  // (gamma-1)

  //-------------------------------------
  // use wrappers to update Qo in-place
  //-------------------------------------
  DVec rho,rhou,rhov,Ener; int Nr=Np*K;
  rho.borrow (Nr, Qo.pCol(1));
  rhou.borrow(Nr, Qo.pCol(2));
  rhov.borrow(Nr, Qo.pCol(3));
  Ener.borrow(Nr, Qo.pCol(4));

  if (1) 
  {
    this->pref = 12.0;

    rho  =  1.0;
    rhou = -sin(2.0*pi*y);
    rhov =  sin(4.0*pi*x);
    Ener = pref/gm1 + 0.5*(sqr(rhou) + sqr(rhov)).dd(rho);
  } 
  else 
  {
    this->pref = 12.0;

    rho  = 1.0;
    rhou = 0.0;
    rhov = 0.0;
    Ener = pref/gm1 + 0.5 * rho.dm(exp(-4.0*(sqr(cos(pi*x))+sqr(cos(pi*y)))));
  }
}

void transpose(const DMat<RT>& A, DMat<RT>& result)
{
  assert(&A != &result);  // these cannot be aliased!
  assert(result.numRows() == A.numColumns());
  assert(result.numColumns() == A.numRows());
  for (size_t c = 0; c < A.numColumns(); ++c)
    {
      auto i = A.columnBegin(c);
      auto j = result.rowBegin(c);
      auto end = A.columnEnd(c);
      for (; i != end; ++i, ++j) A.ring().set(*j, *i);
    }
}

engine_RawArrayIntPairOrNull rawLQUPFactorizationInPlace(MutableMatrix *A, M2_bool transpose)
{
#ifdef HAVE_FFLAS_FFPACK
  // Suppose A is m x n
  // then we get A = LQUP = LSP, see e.g. http://www.ens-lyon.fr/LIP/Pub/Rapports/RR/RR2006/RR2006-28.pdf
  // P and Q are permutation info using LAPACK's convention:, see
  // http://www.netlib.org/lapack/explore-html/d0/d39/_v_a_r_i_a_n_t_s_2lu_2_r_e_c_2dgetrf_8f.html
  // P is n element permutation on column: size(P)=min(m,n); 
  // for 1 <= i <= min(m,n), col i of the  matrix was interchanged with col P(i).
  // Qt is m element permutation on rows (inverse permutation)
  // for 1 <= i <= min(m,n), col i of the  matrix was interchanged with col P(i).
  A->transpose();

  DMat<M2::ARingZZpFFPACK> *mat = A->coerce< DMat<M2::ARingZZpFFPACK> >();
  if (mat == 0) 
    {
      throw exc::engine_error("LUDivine not defined for this ring");
      //      ERROR("LUDivine not defined for this ring");
      //      return 0;
    }
  size_t nelems = mat->numColumns();
  if (mat->numRows() < mat->numColumns()) nelems = mat->numRows();
 

  std::vector<size_t> P(nelems, -1);
  std::vector<size_t> Qt(nelems, -1);
  
  // ignore return value (rank) of:
  LUdivine(mat->ring().field(),
                       FFLAS::FflasNonUnit,
                       (transpose ? FFLAS::FflasTrans : FFLAS::FflasNoTrans),
                       mat->numRows(),
                       mat->numColumns(),
                       mat->array(),
                       mat->numColumns(),
                       &P[0], 
                       &Qt[0]);

  engine_RawArrayIntPairOrNull result = new engine_RawArrayIntPair_struct;
  result->a = stdvector_to_M2_arrayint(Qt);
  result->b = stdvector_to_M2_arrayint(P);

  return result;
#endif
  return 0;
}

//---------------------------------------------------------
DMat& NDG2D::CurvedDGJump2D
(
  const DMat& gU,     // [in]
  const IVec& gmapD,  // [in]
  const DVec& bcU     // [in]
)
//---------------------------------------------------------
{
  // function [jumpU] = CurvedDGJump2D(gU, gmapD, bcU)
  // purpose: compute discontinuous Galerkin jump applied
  //          to a field given at cubature and Gauss nodes

  DMat mmCHOL;  DVec gUM,gUP,fx;
  DMat *tmp = new DMat("jumpU", OBJ_temp);
  DMat &jumpU(*tmp);

  // shorthand references
  Cub2D& cub = this->m_cub; Gauss2D& gauss = this->m_gauss;

  // surface traces at Gauss nodes
  gUM = gU(gauss.mapM);
  gUP = gU(gauss.mapP);
  if (gmapD.size()>0) { gUP(gmapD) = bcU(gmapD); }

  // compute jump term and lift to triangle interiors
  fx = gUM - gUP;
  jumpU = -gauss.interpT*(gauss.W.dm(fx));

  // multiply straight sided triangles by inverse mass matrix
  jumpU(All,straight) = VVT * dd(jumpU(All,straight), J(All,straight));

  // multiply by custom inverse mass matrix for each curvilinear triangle
  int Ncurved = curved.size(), k=0;
  for (int m=1; m<=Ncurved; ++m) {
    k = curved(m);
    mmCHOL.borrow(Np,Np, cub.mmCHOL.pCol(k));
    mmCHOL.set_factmode(FACT_CHOL);  // indicate factored state
    jumpU(All,k) = chol_solve(mmCHOL, jumpU(All,k));
  }

  // these parameters may be OBJ_temp (allocated on the fly)
  if (OBJ_temp == gU.get_mode())  { delete (&gU); }
  if (OBJ_temp == bcU.get_mode()) { delete (&bcU); }

  return jumpU;
}