C++ DoubleMatrix::isEmpty方法代码示例

void multiply(const DoubleMatrix &X, const DoubleMatrix &Y, DoubleMatrix& Z)
{
  assert(&X!=&Z);
  assert(&Y!=&Z);
  if( X.isEmpty() || Y.isEmpty() )
    return;

  //
  // Row Major Matrices
  //
  //    C = alpha * A * B + beta * C  --- (1)
  //
  // A : m by k         lda (stride) = k
  // B : k by n         ldb (stride) = n
  // C : m by n         ldc (stride) = n
  //
  // Column Major Matrices
  //
  //    Z = alpha * X * Y + beta * C  --- (2)
  //    Z = C^t 
  //      = alpha * B^t * A^t + beta * C^t  --- (3)
  //
  // X = B^t : n by k   ldx (stride) = n
  // Y = A^t : k by m   ldy (stride) = k
  // Z = C^t : n by m   ldz (stride) = n
  //
  int m = Y.nc();
  int k = X.nc();
  int n = X.nr();

  Z.resize( n, m );
  
  assert( X.nr() == n );
  assert( X.nc() == k );
  assert( Y.nr() == k );
  assert( Y.nc() == m );
  assert( Z.nr() == n );
  assert( Z.nc() == m );
  
  const double * pX = X.data();
  const double * pY = Y.data();
  double       * pZ = Z.data(); 
  
  int lda = n;
  int ldb = k;
  int ldc = n;

  cblas_dgemm( CblasColMajor, CblasNoTrans, CblasNoTrans, n, m, k, ALPHA, pX, lda, pY, ldb , BETA, pZ, ldc );

}

const DoubleMatrix backDiv(const DoubleMatrix &dmatA, const DoubleMatrix &dmatB)
{
  // A is assumed to be square.
  int m = dmatA.nr();
  int n = dmatA.nc();
  assert( m == n );

  // B is m by l matrix, where l is the number of right hand sides.
  int l = dmatB.nc();
  assert( dmatB.nr() == m );

  if( dmatA.isEmpty() || dmatB.isEmpty() )
    return DoubleMatrix( 0, 0 );

  //==============================================================
  // First decompose A into LU such that A = P * L * U, 
  // where P is the permutation matrix,
  // L is the lower triangle and the U the upper triangle.
  //
  // We use CLAPACK's DGETRF() which does LU decomposition
  // with partial (ie. row interchanges only) pivoting.
  //==============================================================
  
  // enum CBLAS_ORDER order =: (CblasColMajor | CblasRowMajor)
  //
  // If order = CblasColMajor, the array, a, is assumed to 
  // hold each matrix A's column in the contiguous manner
  // in memory (ie. A is said to be in the column major order).
  // If order = CblasRowMajor, the array, a, is assumed to
  // hold each matrix A's row in the contiguous manner
  // in memory (ie. A is said to be in the row major order).
  enum CBLAS_ORDER order = CblasColMajor;

  // double *a
  //
  // (on entry) a points to the elements of matrix A(m,n) 
  // in the column major order if "order" = CblasColMajor, 
  // or in the row major order if "order" = CblasRowMajor.
  //
  // (on exit) The lower triangle (j<=i) is replaced by L
  // and the upper triangle (j>i) is replaced by U.
  double a[m*n];
  copy( dmatA.data(), dmatA.data()+m*n, a );

  // int lda
  // 
  // The leading dimension of A.  
  // If A is in the column major order, lda = m.
  // If A is in the row major order, lda = n.
  int lda = m;

  // int ipiv(m)
  //
  // (on exit) The i-th row in A was interchanged with the row
  // indicated by the value in ipiv[i].
  int ipiv[m];

  int info = clapack_dgetrf( order, m, n, a, lda, ipiv );
  if( info < 0 )
    {
      char mess[ SpkError::maxMessageLen() ];
      snprintf( mess, SpkError::maxMessageLen(), "Solution of a system of linear equations using the LU decomposition failed: \n the %s argument to the function that performs the LU decomposition  had an illegal value.", 
               intToOrdinalString( -info, ONE_IS_FIRST_INT ).c_str() );
      throw SpkException( SpkError::SPK_UNKNOWN_ERR, mess, __LINE__, __FILE__ );
    }
  else if( info > 0 )
    {
      char mess[ SpkError::maxMessageLen() ];
      snprintf( mess, SpkError::maxMessageLen(), "Solution of a system of linear equations using the LU decomposition failed: \nthe %s diagonal element of U is exactly zero.", 
               intToOrdinalString( info, ONE_IS_FIRST_INT ).c_str() );
      throw SpkException( SpkError::SPK_NOT_POS_DEF_ERR, mess, __LINE__, __FILE__ );
    }

  //==============================================================
  // Solve A x = B for x using the LU computed in the previous
  // step.
  // Note that A is now assumed to be square: m = n.
  //==============================================================
  
  // int rhs
  //
  // The number of right hand sides (ie. the number of columns of B).
  int nrhs = l;

  // int ldb
  // The leading dimension of B.
  // If B is in the column major order, ldb = m.
  // If B is in the row major order, ldb = l.
  int ldb = m;

  // double *x
  //
  // (on entry) x points to the elements of B in the column major
  // order if "order" = CblasColMajor or in the row major otherwise.
  // (on exit) x points to the solution matrix, x (m=n by l).
  DoubleMatrix X( dmatB );
  double * x = X.data();// This points to b on entry and contains the solution x upon exit

  info = clapack_dgetrs( order, CblasNoTrans, n, nrhs, a, lda, ipiv, x, ldb );
  if( info < 0 )
//.........这里部分代码省略.........