C++ SymmetricMatrix类代码示例

void test1(Real* y, Real* x1, Real* x2, int nobs, int npred)
{
   cout << "\n\nTest 1 - traditional, bad\n";

   // traditional sum of squares and products method of calculation
   // but not adjusting means; maybe subject to round-off error

   // make matrix of predictor values with 1s into col 1 of matrix
   int npred1 = npred+1;        // number of cols including col of ones.
   Matrix X(nobs,npred1);
   X.Column(1) = 1.0;

   // load x1 and x2 into X
   //    [use << rather than = when loading arrays]
   X.Column(2) << x1;  X.Column(3) << x2;

   // vector of Y values
   ColumnVector Y(nobs); Y << y;

   // form sum of squares and product matrix
   //    [use << rather than = for copying Matrix into SymmetricMatrix]
   SymmetricMatrix SSQ; SSQ << X.t() * X;

   // calculate estimate
   //    [bracket last two terms to force this multiplication first]
   //    [ .i() means inverse, but inverse is not explicity calculated]
   ColumnVector A = SSQ.i() * (X.t() * Y);

   // Get variances of estimates from diagonal elements of inverse of SSQ
   // get inverse of SSQ - we need it for finding D
   DiagonalMatrix D; D << SSQ.i();
   ColumnVector V = D.AsColumn();

   // Calculate fitted values and residuals
   ColumnVector Fitted = X * A;
   ColumnVector Residual = Y - Fitted;
   Real ResVar = Residual.SumSquare() / (nobs-npred1);

   // Get diagonals of Hat matrix (an expensive way of doing this)
   DiagonalMatrix Hat;  Hat << X * (X.t() * X).i() * X.t();

   // print out answers
   cout << "\nEstimates and their standard errors\n\n";

   // make vector of standard errors
   ColumnVector SE(npred1);
   for (int i=1; i<=npred1; i++) SE(i) = sqrt(V(i)*ResVar);
   // use concatenation function to form matrix and use matrix print
   // to get two columns
   cout << setw(11) << setprecision(5) << (A | SE) << endl;

   cout << "\nObservations, fitted value, residual value, hat value\n";

   // use concatenation again; select only columns 2 to 3 of X
   cout << setw(9) << setprecision(3) <<
     (X.Columns(2,3) | Y | Fitted | Residual | Hat.AsColumn());
   cout << "\n\n";
}

//--------------------------------------------
void NRLib::CholeskyInvert(SymmetricMatrix& A)
//--------------------------------------------
{
  //NBNB legge til feilmld her, se metode ovenfor
  int infof = flens::potrf(A.upLo(), A.dim(), A.data(), A.leadingDimension());
  int infoi = flens::potri(A.upLo(), A.dim(), A.data(), A.leadingDimension());

  if(infof > 0 || infoi > 0)
    throw NRLib::Exception("Error in Cholesky inversion");
}

  inline
  Matrix<C>::Matrix(const SymmetricMatrix<C>& M)
    : entries_(M.row_dimension()*M.column_dimension()),
      rowdim_(M.row_dimension()), coldim_(M.column_dimension())
  {
    for (size_type i(0); i < rowdim_; i++)
      for (size_type j(0); j <= i; j++)
	{
	  this->operator () (i, j) = this->operator () (j, i) = M(i, j);
	}
  }

void SpectClust::normalizeSum(SymmetricMatrix &A) {
  RowVector Ones(A.Ncols());
  Ones = 1.0;
  // Calculate the sum of each row matrix style. Could this be a Diagonal Matrix instead?
  Matrix RowSum = Ones * (A.t());
  DiagonalMatrix D(A.Ncols());
  D = 0;
  // Take the inverse square root
  for(int i = 1; i <= A.Ncols(); i++) {
    D(i,i) = 1/sqrt(RowSum(1,i));
  }
  Matrix X = D * (A * D);
  A << X;
}

std::vector<int> fiedler_reorder(const SymmetricMatrix& m)
{
  SymmetricMatrix absm=m;
  const int nrows=m.Nrows();
  for (int i=0;i<nrows;++i) {
    for (int j=0;j<=i;++j){
       //absolute value
       absm.element(i,j)=std::fabs(absm.element(i,j));
    }
  }

  //laplacian
  SymmetricMatrix lap(nrows);
  lap=0.;
  for (int i=0;i<nrows;++i)
    lap.element(i,i)=absm.Row(i+1).Sum();
  lap-=absm;

  DiagonalMatrix eigs; 
  Matrix vecs;
  
  EigenValues(lap,eigs,vecs);

  ColumnVector fvec=vecs.Column(2);
  std::vector<double> fvec_stl(nrows);
  //copies over fvec to fvec_stl
  std::copy(&fvec.element(0),&fvec.element(0)+nrows,fvec_stl.begin());
  std::vector<int> findices;
  //sorts the data by eigenvalue in ascending order
  sort_data_to_indices(fvec_stl,findices);
  
  return findices;
  /* BLOCK works with findices*/

}

void SpectClust::fillInDistance(SymmetricMatrix &A, const Matrix &M, const DistanceMetric &dMetric, bool expon) {
  A.ReSize(M.Ncols());
  for(int i = 1; i <= M.Ncols(); i++) {
    for(int j = i; j <= M.Ncols(); j++) {
      if(expon) 
        A(j,i) = A(i,j) = exp(-1 * dMetric.dist(M,i,j));
      else
        A(j,i) = A(i,j) = dMetric.dist(M,i,j);
    }
  }
}

static void tred3(const SymmetricMatrix& X, DiagonalMatrix& D,
   DiagonalMatrix& E, SymmetricMatrix& A)
{
   Tracer et("Evalue(tred3)");
   Real tol =
      FloatingPointPrecision::Minimum()/FloatingPointPrecision::Epsilon();
   int n = X.Nrows(); A = X; D.ReSize(n); E.ReSize(n);
   Real* ei = E.Store() + n;
   for (int i = n-1; i >= 0; i--)
   {
      Real h = 0.0; Real f;
      Real* d = D.Store(); Real* a = A.Store() + (i*(i+1))/2; int k = i;
      while (k--) { f = *a++; *d++ = f; h += square(f); }
      if (h <= tol) { *(--ei) = 0.0; h = 0.0; }
      else
      {
	 Real g = sign(-sqrt(h), f); *(--ei) = g; h -= f*g;
         f -= g; *(d-1) = f; *(a-1) = f; f = 0.0;
         Real* dj = D.Store(); Real* ej = E.Store(); int j;
         for (j = 0; j < i; j++)
         {
            Real* dk = D.Store(); Real* ak = A.Store()+(j*(j+1))/2;
            Real g = 0.0; k = j;
            while (k--)  g += *ak++ * *dk++;
            k = i-j; int l = j; 
            while (k--) { g += *ak * *dk++; ak += ++l; }
	    g /= h; *ej++ = g; f += g * *dj++;
         }  
	 Real hh = f / (2 * h); Real* ak = A.Store();
         dj = D.Store(); ej = E.Store();
         for (j = 0; j < i; j++)
         {
	    f = *dj++; g = *ej - hh * f; *ej++ = g;
            Real* dk = D.Store(); Real* ek = E.Store(); k = j+1;
	    while (k--) { *ak++ -= (f * *ek++ + g * *dk++); }
	 }
      }
      *d = *a; *a = h;
   }
}

ReturnMatrix Cholesky(const SymmetricMatrix& S)
{
   REPORT
   Tracer trace("Cholesky");
   int nr = S.Nrows();
   LowerTriangularMatrix T(nr);
   Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
   for (int i=0; i<nr; i++)
   {
      Real* tj = t; Real sum; int k;
      for (int j=0; j<i; j++)
      {
         Real* tk = ti; sum = 0.0; k = j;
         while (k--) { sum += *tj++ * *tk++; }
         *tk = (*s++ - sum) / *tj++;
      }
      sum = 0.0; k = i;
      while (k--) { sum += square(*ti++); }
      Real d = *s++ - sum;
      if (d<=0.0)  Throw(NPDException(S));
      *ti++ = sqrt(d);
   }
   T.release(); return T.for_return();
}

bool SpectClust::findNLargestSymEvals(const SymmetricMatrix &W, int numLamda, std::vector<Numeric> &eVals, Matrix &EVec) {
  bool converged = false;
  eVals.clear();
  eVals.reserve(numLamda);
  DiagonalMatrix D(W.Ncols());
  Matrix E;
  try {
    EigenValues(W, D, E);
    converged = true;
    EVec.ReSize(W.Ncols(), numLamda);
    int count = 0;
    for(int i = W.Ncols(); i > W.Ncols() - numLamda; i--) {
      eVals.push_back(D(i));
      EVec.Column(++count) << E.Column(i);
    }
  }
  catch(const Exception &e) {
    Err::errAbort("Exception: " + ToStr(e.what()));
  }
  catch(...) {
    Err::errAbort("Yikes couldn't calculate eigen vectors.");
  }
  return converged;
}

int
AccpmLASymmLinSolve(SymmetricMatrix &A, RealMatrix &X, const RealMatrix &B)
{
  LaLowerTriangMatDouble L(A);
  char uplo = 'L';
  long int n = A.size(0);
  long int lda = A.inc(0) * A.gdim(0);
  long int nrhs = B.size(1);
  long int ldb = B.inc(0) * B.gdim(0);
  long int info;
  
  X.inject(B);

  F77NAME(dposv) (&uplo, &n, &nrhs, &L(0,0), &lda, &X(0,0), &ldb, &info);
 
  if (info > 0) {
    std::cerr << "AccpmLASymmLinSolve: Symmetric Matrix is not positive-definite." 
	      << std::endl;
  } else if (info < 0) {
    std::cerr << "AccpmLASymmLinSolve: argument " << -info << " is invalid" 
	      << std::endl;
  }
  return info;
}

void Print(const SymmetricMatrix& X)
{
   ++PCN;
   cout << "\nMatrix type: " << X.Type().Value() << " (";
   cout << X.Nrows() << ", ";
   cout << X.Ncols() << ")\n\n";
   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
   int nr=X.Nrows(); int nc=X.Ncols();
   for (int i=1; i<=nr; i++)
   {
      int j;
      for (j=1; j<i; j++) cout << X(j,i) << "\t";
      for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
      cout << "\n";
   }
   cout << flush; ++PCZ;
}

void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
   Matrix& V, bool eivec)
{
   Real epsilon = FloatingPointPrecision::Epsilon();
   Tracer et("Jacobi");
   REPORT
   int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.resize(n); A = X;
   if (eivec) { REPORT V.resize(n,n); D = 1.0; V = D; }
   B << A; D = B; Z = 0.0; A.Inject(Z);
   bool converged = false;
   for (int i=1; i<=50; i++)
   {
      Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
      while (p--) sm += fabs(*a++);            // have previously zeroed diags
      if (sm==0.0) { REPORT converged = true; break; }
      Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
      for (p = 0; p < n; p++)
      {
         Real* ap1 = a + (p*(p+1))/2;
         Real& zp = Z.element(p); Real& dp = D.element(p);
         for (int q = p+1; q < n; q++)
         {
            Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
            Real& zq = Z.element(q); Real& dq = D.element(q);
            Real& apq = A.element(q,p);
            Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);

            if (i>4 && g < epsilon*adp && g < epsilon*adq) { REPORT apq = 0.0; }
            else if (fabs(apq) > tresh)
            {
               REPORT
               Real t; Real h = dq - dp; Real ah = fabs(h);
               if (g < epsilon*ah) { REPORT t = apq / h; }
               else
               {
                  REPORT
                  Real theta = 0.5 * h / apq;
                  t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
                  if (theta<0.0) { REPORT t = -t; }
               }
               Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
               Real tau = s / (1.0 + c); h = t * apq;
               zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
               int j = p;
               while (j--)
               {
                  g = *ap; h = *aq;
                  *ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
               }
               int ip = p+1; j = q-ip; ap += ip++; aq++;
               while (j--)
               {
                  g = *ap; h = *aq;
                  *ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
                  ap += ip++;
               }
               if (q < n-1)             // last loop is non-empty
               {
                  int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
                  for (;;)
                  {
                     g = *ap; h = *aq;
                     *ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
                     if (!(--j)) break;
                     ap += ip++; aq += iq++;
                  }
               }
               if (eivec)
               {
                  REPORT
                  RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
                  Rotate(VP, VQ, tau, s);
               }
            }
         }
      }
      B = B + Z; D = B; Z = 0.0;
   }
   if (!converged) Throw(ConvergenceException(X));
   if (eivec) SortSV(D, V, true);
   else SortAscending(D);
}

/*
 *	s o l v e O q p B e n c h m a r k
 */
returnValue solveOqpBenchmark(	int_t nQP, int_t nV,
								const real_t* const _H, const real_t* const g,
								const real_t* const lb, const real_t* const ub,
								BooleanType isSparse, BooleanType useHotstarts, 
								const Options& options, int_t maxAllowedNWSR,
								real_t& maxNWSR, real_t& avgNWSR, real_t& maxCPUtime, real_t& avgCPUtime,
								real_t& maxStationarity, real_t& maxFeasibility, real_t& maxComplementarity
								)
{
	int_t k;

	/* I) SETUP AUXILIARY VARIABLES: */
	/* 1) Keep nWSR and store current and maximum number of
	 *    working set recalculations in temporary variables */
	int_t nWSRcur;

	real_t CPUtimeLimit = maxCPUtime;
	real_t CPUtimeCur = CPUtimeLimit;
	real_t stat, feas, cmpl;
	maxNWSR = 0;
	avgNWSR = 0;
	maxCPUtime = 0.0;
	avgCPUtime = 0.0;
	maxStationarity = 0.0;
	maxFeasibility = 0.0;
	maxComplementarity = 0.0;

	/* 2) Pointers to data of current QP ... */
	const real_t* gCur;
	const real_t* lbCur;
	const real_t* ubCur;

	/* 3) Vectors for solution obtained by qpOASES. */
	real_t* x = new real_t[nV];
	real_t* y = new real_t[nV];
	//real_t  obj;

	/* 4) Prepare matrix objects */
	SymmetricMatrix *H; 
	real_t* H_cpy = new real_t[nV*nV];
	memcpy( H_cpy,_H, ((uint_t)(nV*nV))*sizeof(real_t) );

	if ( isSparse == BT_TRUE )
	{
		SymSparseMat *Hs;
		H = Hs = new SymSparseMat(nV, nV, nV, H_cpy);
		Hs->createDiagInfo();
		delete[] H_cpy;
	}
	else
	{
		H = new SymDenseMat(nV, nV, nV, const_cast<real_t *>(H_cpy));
	}
	
	H->doFreeMemory( );

	/* II) SETUP QPROBLEM OBJECT */
	QProblemB qp( nV );
	qp.setOptions( options );
	//qp.setPrintLevel( PL_LOW );


	/* III) RUN BENCHMARK SEQUENCE: */
	returnValue returnvalue;

	for( k=0; k<nQP; ++k )
	{
		//if ( k%50 == 0 )
		//	printf( "%d\n",k );

		/* 1) Update pointers to current QP data. */
		gCur   = &( g[k*nV] );
		lbCur  = &( lb[k*nV] );
		ubCur  = &( ub[k*nV] );

		/* 2) Set nWSR and maximum CPU time. */
		nWSRcur = maxAllowedNWSR;
		CPUtimeCur = CPUtimeLimit;

		/* 3) Solve current QP. */
		if ( ( k == 0 ) || ( useHotstarts == BT_FALSE ) )
		{
			/* initialise */
			returnvalue = qp.init( H,gCur,lbCur,ubCur, nWSRcur,&CPUtimeCur );
			if ( ( returnvalue != SUCCESSFUL_RETURN ) && ( returnvalue != RET_MAX_NWSR_REACHED ) )
			{
				delete H; delete[] y; delete[] x;
				return THROWERROR( returnvalue );
			}
		}
		else
		{
			/* hotstart */
			returnvalue = qp.hotstart( gCur,lbCur,ubCur, nWSRcur,&CPUtimeCur );
			if ( ( returnvalue != SUCCESSFUL_RETURN ) && ( returnvalue != RET_MAX_NWSR_REACHED ) )
			{
				delete H; delete[] y; delete[] x;
//.........这里部分代码省略.........

void trymatd()
{
   Tracer et("Thirteenth test of Matrix package");
   Tracer::PrintTrace();
   Matrix X(5,20);
   int i,j;
   for (j=1;j<=20;j++) X(1,j) = j+1;
   for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
   SymmetricMatrix S; S << X * X.t();
   Matrix SM = X * X.t() - S;
   Print(SM);
   LowerTriangularMatrix L = Cholesky(S);
   Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
   Print(Diff);
   {
      Tracer et1("Stage 1");
      LowerTriangularMatrix L1(5);
      Matrix Xt = X.t(); Matrix Xt2 = Xt;
      QRZT(X,L1);
      Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
      UpperTriangularMatrix Ut(5);
      QRZ(Xt,Ut);
      Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
      Matrix Y(3,20);
      for (j=1;j<=20;j++) Y(1,j) = 22-j;
      for (i=2;i<=3;i++) for (j=1;j<=20; j++)
         Y(i,j) = (long)Y(i-1,j) * j % 101;
      Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
      QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
      Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Y2 * Xt2 * S.i() - M * L.i();
      Clean(Diff,0.000000001); Print(Diff);
   }

   ColumnVector C1(5);
   {
      Tracer et1("Stage 2");
      X.ReSize(5,5);
      for (j=1;j<=5;j++) X(1,j) = j+1;
      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      for (i=1;i<=5;i++) C1(i) = i*i;
      CroutMatrix A = X;
      ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
      X = C1 - C2; Clean(X,0.000000001); Print(X);
   }

   {
      Tracer et1("Stage 3");
      X.ReSize(7,7);
      for (j=1;j<=7;j++) X(1,j) = j+1;
      for (i=2;i<=7;i++) for (j=1;j<=7; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      C1.ReSize(7);
      for (i=1;i<=7;i++) C1(i) = i*i;
      RowVector R1 = C1.t();
      Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
      Print(Diff);
   }

   {
      Tracer et1("Stage 4");
      X.ReSize(5,5);
      for (j=1;j<=5;j++) X(1,j) = j+1;
      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      C1.ReSize(5);
      for (i=1;i<=5;i++) C1(i) = i*i;
      CroutMatrix A1 = X*X;
      ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
      X = C1 - C2; Clean(X,0.000000001); Print(X);
   }


   {
      Tracer et1("Stage 5");
      int n = 40;
      SymmetricBandMatrix B(n,2); B = 0.0;
      for (i=1; i<=n; i++)
      {
         B(i,i) = 6;
         if (i<=n-1) B(i,i+1) = -4;
         if (i<=n-2) B(i,i+2) = 1;
      }
      B(1,1) = 5; B(n,n) = 5;
      SymmetricMatrix A = B;
      ColumnVector X(n);
      X(1) = 429;
      for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
      X = X / 100000L;
      // the matrix B is rather ill-conditioned so the difficulty is getting
      // good agreement (we have chosen X very small) may not be surprising;
      // maximum element size in B.i() is around 1400
      ColumnVector Y1 = A.i() * X;
      LowerTriangularMatrix C1 = Cholesky(A);
      ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
      Clean(Y2, 0.000000001); Print(Y2);
      UpperTriangularMatrix CU = C1.t().i();
//.........这里部分代码省略.........

void trymatd()
{
   Tracer et("Thirteenth test of Matrix package");
   Tracer::PrintTrace();
   Matrix X(5,20);
   int i,j;
   for (j=1;j<=20;j++) X(1,j) = j+1;
   for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
   SymmetricMatrix S; S << X * X.t();
   Matrix SM = X * X.t() - S;
   Print(SM);
   LowerTriangularMatrix L = Cholesky(S);
   Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
   Print(Diff);
   {
      Tracer et1("Stage 1");
      LowerTriangularMatrix L1(5);
      Matrix Xt = X.t(); Matrix Xt2 = Xt;
      QRZT(X,L1);
      Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
      UpperTriangularMatrix Ut(5);
      QRZ(Xt,Ut);
      Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
      Matrix Y(3,20);
      for (j=1;j<=20;j++) Y(1,j) = 22-j;
      for (i=2;i<=3;i++) for (j=1;j<=20; j++)
         Y(i,j) = (long)Y(i-1,j) * j % 101;
      Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
      QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
      Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
      Diff = Y2 * Xt2 * S.i() - M * L.i();
      Clean(Diff,0.000000001); Print(Diff);
   }

   ColumnVector C1(5);
   {
      Tracer et1("Stage 2");
      X.ReSize(5,5);
      for (j=1;j<=5;j++) X(1,j) = j+1;
      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      for (i=1;i<=5;i++) C1(i) = i*i;
      CroutMatrix A = X;
      ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
      X = C1 - C2; Clean(X,0.000000001); Print(X);
   }

   {
      Tracer et1("Stage 3");
      X.ReSize(7,7);
      for (j=1;j<=7;j++) X(1,j) = j+1;
      for (i=2;i<=7;i++) for (j=1;j<=7; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      C1.ReSize(7);
      for (i=1;i<=7;i++) C1(i) = i*i;
      RowVector R1 = C1.t();
      Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
      Print(Diff);
   }

   {
      Tracer et1("Stage 4");
      X.ReSize(5,5);
      for (j=1;j<=5;j++) X(1,j) = j+1;
      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
         X(i,j) = (long)X(i-1,j) * j % 1001;
      C1.ReSize(5);
      for (i=1;i<=5;i++) C1(i) = i*i;
      CroutMatrix A1 = X*X;
      ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
      X = C1 - C2; Clean(X,0.000000001); Print(X);
   }


   {
      Tracer et1("Stage 5");
      int n = 40;
      SymmetricBandMatrix B(n,2); B = 0.0;
      for (i=1; i<=n; i++)
      {
         B(i,i) = 6;
         if (i<=n-1) B(i,i+1) = -4;
         if (i<=n-2) B(i,i+2) = 1;
      }
      B(1,1) = 5; B(n,n) = 5;
      SymmetricMatrix A = B;
      ColumnVector X(n);
      X(1) = 429;
      for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
      X = X / 100000L;
      // the matrix B is rather ill-conditioned so the difficulty is getting
      // good agreement (we have chosen X very small) may not be surprising;
      // maximum element size in B.i() is around 1400
      ColumnVector Y1 = A.i() * X;
      LowerTriangularMatrix C1 = Cholesky(A);
      ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
      Clean(Y2, 0.000000001); Print(Y2);
      UpperTriangularMatrix CU = C1.t().i();
//.........这里部分代码省略.........