C++ SymmetricMatrix::ReSize方法代码示例

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

void trymat7()
{
//   cout << "\nSeventh test of Matrix package\n";
    Tracer et("Seventh test of Matrix package");
    Tracer::PrintTrace();

    int i,j;


    DiagonalMatrix D(6);
    UpperTriangularMatrix U(6);
    for (i=1; i<=6; i++) {
        for (j=i; j<=6; j++) U(i,j)=i*i*j-50;
        D(i,i)=i*i+i-10;
    }
    LowerTriangularMatrix L=(U*3.0).t();
    SymmetricMatrix S(6);
    for (i=1; i<=6; i++) for (j=i; j<=6; j++) S(i,j)=i*i+2.0+j;
    Matrix MD=D;
    Matrix ML=L;
    Matrix MU=U;
    Matrix MS=S;
    Matrix M(6,6);
    for (i=1; i<=6; i++) for (j=1; j<=6; j++) M(i,j)=i*j+i*i-10.0;
    {
        Tracer et1("Stage 1");
        Print(Matrix((S-M)-(MS-M)));
        Print(Matrix((-M-S)+(MS+M)));
        Print(Matrix((U-M)-(MU-M)));
    }
    {
        Tracer et1("Stage 2");
        Print(Matrix((L-M)+(M-ML)));
        Print(Matrix((D-M)+(M-MD)));
        Print(Matrix((D-S)+(MS-MD)));
        Print(Matrix((D-L)+(ML-MD)));
    }

    {
        M=MU.t();
    }
    LowerTriangularMatrix LY=D.i()*U.t();
    {
        Tracer et1("Stage 3");
        MS=D*LY-M;
        Clean(MS,0.00000001);
        Print(MS);
        L=U.t();
        LY=D.i()*L;
        MS=D*LY-M;
        Clean(MS,0.00000001);
        Print(MS);
    }
    {
        Tracer et1("Stage 4");
        UpperTriangularMatrix UT(11);
        int i, j;
        for (i=1; i<=11; i++) for (j=i; j<=11; j++) UT(i,j)=i*i+j*3;
        GenericMatrix GM;
        Matrix X;
        UpperBandMatrix UB(11,3);
        UB.Inject(UT);
        UT = UB;
        UpperBandMatrix UB2 = UB / 8;
        GM = UB2-UT/8;
        X = GM;
        Print(X);
        SymmetricBandMatrix SB(11,4);
        SB << (UB + UB.t());
        X = SB - UT - UT.t();
        Print(X);
        BandMatrix B = UB + UB.t()*2;
        DiagonalMatrix D;
        D << B;
        X.ReSize(1,1);
        X(1,1) = Trace(B)-Sum(D);
        Print(X);
        X = SB + 5;
        Matrix X1=X;
        X = SP(UB,X);
        Matrix X2 =UB;
        X1 = (X1.AsDiagonal() * X2.AsDiagonal()).AsRow()-X.AsColumn().t();
        Print(X1);
        X1=SB.t();
        X2 = B.t();
        X = SB.i() * B - X1.i() * X2.t();
        Clean(X,0.00000001);
        Print(X);
        X = SB.i();
        X = X * B - X1.i() * X2.t();
        Clean(X,0.00000001);
        Print(X);
        D = 1;
        X = SB.i() * SB - D;
        Clean(X,0.00000001);
        Print(X);
        ColumnVector CV(11);
        CV << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
        D << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
        X = CV.AsDiagonal();
//.........这里部分代码省略.........

void init_illum (int ndim, ColumnVector& x)
{
  int    i, j, k;
  double h,rij2,midx,midy,dx,dy,slope1,slope2,theta1,theta2,theta,scale;
  double dtmp, Ides;

  // allocate storage

  m = ndim;
  n = 11;
  Ides = 2.0;
  patch = new double*[n+1];
  for (i=0; i<=n; i++) patch[i] = new double[2];
  lamp = new double*[m+1];
  for (i=0; i<=m; i++) lamp[i] = new double[2];
  A = new double*[n+1];
  for (i=0; i<=n; i++) A[i] = new double[m+1];

  // initializing the patches and lamps

  patch[0][0] = 0.0;     patch[0][1] = 0.0;
  patch[1][0] = 0.0909;  patch[1][1] = 0.1;
  patch[2][0] = 0.1818;  patch[2][1] = 0.2;
  patch[3][0] = 0.2727;  patch[3][1] = 0.2;
  patch[4][0] = 0.3636;  patch[4][1] = 0.1;
  patch[5][0] = 0.4545;  patch[5][1] = 0.2;
  patch[6][0] = 0.5455;  patch[6][1] = 0.3;
  patch[7][0] = 0.6364;  patch[7][1] = 0.2;
  patch[8][0] = 0.7273;  patch[8][1] = 0.0;
  patch[9][0] = 0.8182;  patch[9][1] = 0.0;
  patch[10][0] = 0.9091; patch[10][1] = 0.2;
  patch[11][0] = 1.0;    patch[11][1] = 0.1;

  lamp[1][0] = 0.1;  lamp[1][1] = 1.0;
  lamp[2][0] = 0.3;  lamp[2][1] = 1.1;
  lamp[3][0] = 0.4;  lamp[3][1] = 0.6;
  lamp[4][0] = 0.6;  lamp[4][1] = 0.9;
  lamp[5][0] = 0.8;  lamp[5][1] = 0.9;
  lamp[6][0] = 0.9;  lamp[6][1] = 1.2;
  lamp[7][0] = 0.95; lamp[7][1] = 1.0;

  // initialize the A matrix

  for (i=1; i<=n; i++) {
    midx  = patch[i][0] - patch[i-1][0];
    midy  = patch[i][1] - patch[i-1][1];
    if (midx == 0.0)  {
      cout << "Error : patches cannot overlap each other. \n";
      exit(-1);
    }
    slope1 = midy / midx;
    theta1 = atan(slope1);
    midx  = 0.5 * midx + patch[i-1][0];
    midy  = 0.5 * midy + patch[i-1][1];
    for (j=1; j<=m; j++) {
      dx = lamp[j][0] - midx;
      dy = lamp[j][1] - midy;
      rij2 = dx * dx + dy * dy;
      if (dx == 0.0) theta2 = dy / fabs(dy) * M_PI * 0.5;
      else {
        slope2 = dy / dx;
        theta2 = atan(slope2);
      }
      if (theta2 < 0.0) theta2 = theta2 + M_PI;
      theta = M_PI * 0.5 - (theta2 - theta1);
      scale = cos(theta);
      if (scale < 0.0) A[i][j] = 0.0;
      else A[i][j] = scale / (rij2 * Ides);
    }
  }

  // initialize x 

  for (i=1; i<=m; i++) x(i) = 0.5;

  // initialize the Hessian

  HH.ReSize(m);
  for (i=1; i<=m; i++) {
    for (j=1; j<=m; j++) {
      dtmp = 0.0;
      for (k=1; k<=n; k++) dtmp += A[k][i] * A[k][j];
      HH(i,j) = dtmp;
    }
  }
}