C++ Matrix函数代码示例

MeshData * JSONDataParser::_parseMesh(const rapidjson::Value & rawData)
{
    const auto mesh = BaseObject::borrowObject<MeshData>();

    const auto rawVertices = rawData[VERTICES].GetArray();
    const auto rawUVs = rawData[UVS].GetArray();
    const auto rawTriangles = rawData[TRIANGLES].GetArray();

    const auto numVertices = (unsigned)(rawVertices.Size() / 2);
    const auto numTriangles = (unsigned)(rawTriangles.Size() / 3);

    std::vector<Matrix> inverseBindPose(this->_armature->getSortedBones().size(), Matrix());

    mesh->skinned = rawData.HasMember(WEIGHTS) && !rawData[WEIGHTS].GetArray().Empty();
    mesh->uvs.resize(numVertices * 2);
    mesh->vertices.resize(numVertices * 2);
    mesh->vertexIndices.resize(numTriangles * 3);

    if (mesh->skinned)
    {
        mesh->boneIndices.resize(numVertices);
        mesh->weights.resize(numVertices);
        mesh->boneVertices.resize(numVertices);

        if (rawData.HasMember(SLOT_POSE))
        {
            const auto& rawSlotPose = rawData[SLOT_POSE].GetArray();
            mesh->slotPose.a = rawSlotPose[0].GetFloat();
            mesh->slotPose.b = rawSlotPose[1].GetFloat();
            mesh->slotPose.c = rawSlotPose[2].GetFloat();
            mesh->slotPose.d = rawSlotPose[3].GetFloat();
            mesh->slotPose.tx = rawSlotPose[4].GetFloat();
            mesh->slotPose.ty = rawSlotPose[5].GetFloat();
        }

        if (rawData.HasMember(BONE_POSE))
        {
            const auto& rawBonePose = rawData[BONE_POSE].GetArray();
            for (std::size_t i = 0, l = rawBonePose.Size(); i < l; i += 7)
            {
                const auto rawBoneIndex = rawBonePose[i].GetUint();
                auto& boneMatrix = inverseBindPose[rawBoneIndex];
                boneMatrix.a = rawBonePose[i + 1].GetFloat();
                boneMatrix.b = rawBonePose[i + 2].GetFloat();
                boneMatrix.c = rawBonePose[i + 3].GetFloat();
                boneMatrix.d = rawBonePose[i + 4].GetFloat();
                boneMatrix.tx = rawBonePose[i + 5].GetFloat();
                boneMatrix.ty = rawBonePose[i + 6].GetFloat();
                boneMatrix.invert();
            }
        }
    }

    for (std::size_t i = 0, iW = 0, l = rawVertices.Size(); i < l; i += 2)
    {
        const auto iN = i + 1;
        const auto vertexIndex = i / 2;

        auto x = mesh->vertices[i] = rawVertices[i].GetFloat() * this->_armatureScale;
        auto y = mesh->vertices[iN] = rawVertices[iN].GetFloat() * this->_armatureScale;
        mesh->uvs[i] = rawUVs[i].GetFloat();
        mesh->uvs[iN] = rawUVs[iN].GetFloat();

        if (mesh->skinned)
        {
            const auto rawWeights = rawData[WEIGHTS].GetArray();
            const auto numBones = rawWeights[iW].GetUint();
            auto& indices = mesh->boneIndices[vertexIndex];
            auto& weights = mesh->weights[vertexIndex];
            auto& boneVertices = mesh->boneVertices[vertexIndex];
            mesh->slotPose.transformPoint(x, y, _helpPoint);
            x = mesh->vertices[i] = _helpPoint.x;
            y = mesh->vertices[iN] = _helpPoint.y;

            for (std::size_t iB = 0; iB < numBones; ++iB)
            {
                const auto iI = iW + 1 + iB * 2;
                const auto rawBoneIndex = rawWeights[iI].GetUint();
                const auto boneData = this->_rawBones[rawBoneIndex];

                const auto iderator = std::find(mesh->bones.cbegin(), mesh->bones.cend(), boneData);
                std::size_t boneIndex = 0;
                if (iderator == mesh->bones.cend())
                {
                    boneIndex = mesh->bones.size();
                    mesh->bones.push_back(boneData);
                    mesh->inverseBindPose.push_back(inverseBindPose[rawBoneIndex]); // copy
                }
                else
                {
                    boneIndex = std::distance(mesh->bones.cbegin(), iderator);
                }

                mesh->inverseBindPose[boneIndex].transformPoint(x, y, _helpPoint);

                indices.push_back(boneIndex);
                weights.push_back(rawWeights[iI + 1].GetFloat());
                boneVertices.push_back(_helpPoint.x);
                boneVertices.push_back(_helpPoint.y);
            }
//.........这里部分代码省略.........

MeshComponent::MeshComponent(Mesh* mesh, float scale) {
	this->mesh = mesh;
	this->scale = Matrix();
	this->scale *= scale;
}

Response*
ElasticForceBeamColumn2d::setResponse(const char **argv, int argc, OPS_Stream &output)
{
  Response *theResponse = 0;

  output.tag("ElementOutput");
  output.attr("eleType","ElasticForceBeamColumn2d");
  output.attr("eleTag",this->getTag());
  output.attr("node1",connectedExternalNodes[0]);
  output.attr("node2",connectedExternalNodes[1]);

  // global force - 
  if (strcmp(argv[0],"forces") == 0 || strcmp(argv[0],"force") == 0
      || strcmp(argv[0],"globalForce") == 0 || strcmp(argv[0],"globalForces") == 0) {

    output.tag("ResponseType","Px_1");
    output.tag("ResponseType","Py_1");
    output.tag("ResponseType","Mz_1");
    output.tag("ResponseType","Px_2");
    output.tag("ResponseType","Py_2");
    output.tag("ResponseType","Mz_2");

    theResponse =  new ElementResponse(this, 1, theVector);
  
  
  // local force -
  } else if (strcmp(argv[0],"localForce") == 0 || strcmp(argv[0],"localForces") == 0) {

    output.tag("ResponseType","N_1");
    output.tag("ResponseType","V_1");
    output.tag("ResponseType","M_1");
    output.tag("ResponseType","N_2");
    output.tag("ResponseType","V_2");
    output.tag("ResponseType","M_2");

    theResponse =  new ElementResponse(this, 2, theVector);
  

  // basic force -
  } else if (strcmp(argv[0],"basicForce") == 0 || strcmp(argv[0],"basicForces") == 0) {

    output.tag("ResponseType","N");
    output.tag("ResponseType","M_1");
    output.tag("ResponseType","M_2");

    theResponse =  new ElementResponse(this, 7, Vector(3));

  // chord rotation -
  } else if (strcmp(argv[0],"chordRotation") == 0 || strcmp(argv[0],"chordDeformation") == 0 
	     || strcmp(argv[0],"basicDeformation") == 0) {

    output.tag("ResponseType","eps");
    output.tag("ResponseType","theta_1");
    output.tag("ResponseType","theta_2");

    theResponse =  new ElementResponse(this, 3, Vector(3));
  
  // plastic rotation -
  } else if (strcmp(argv[0],"plasticRotation") == 0 || strcmp(argv[0],"plasticDeformation") == 0) {

    output.tag("ResponseType","epsP");
    output.tag("ResponseType","thetaP_1");
    output.tag("ResponseType","thetaP_2");

    theResponse =  new ElementResponse(this, 4, Vector(3));

  // point of inflection
  } else if (strcmp(argv[0],"inflectionPoint") == 0) {
    
    output.tag("ResponseType","inflectionPoint");

    theResponse =  new ElementResponse(this, 5, 0.0);
  
  // tangent drift
  } else if (strcmp(argv[0],"tangentDrift") == 0) {
    theResponse =  new ElementResponse(this, 6, Vector(2));

  // basic forces
  } else if (strcmp(argv[0],"basicForce") == 0)
    theResponse = new ElementResponse(this, 7, Vector(3));

  else if (strcmp(argv[0],"integrationPoints") == 0)
    theResponse = new ElementResponse(this, 10, Vector(numSections));

  else if (strcmp(argv[0],"integrationWeights") == 0)
    theResponse = new ElementResponse(this, 11, Vector(numSections));

  else if (strcmp(argv[0],"basicStiffness") == 0)
    theResponse = new ElementResponse(this, 12, Matrix(3,3));

  // section response -
  else if (strstr(argv[0],"sectionX") != 0) {
    if (argc > 2) {
      double sectionLoc = atof(argv[1]);

      double xi[maxNumSections];
      double L = crdTransf->getInitialLength();
      beamIntegr->getSectionLocations(numSections, L, xi);
      
      sectionLoc /= L;
//.........这里部分代码省略.........

Matrix Matrix::operator/(float f) const
{
    return Matrix(m_Matrix / f);
}

int main (int argc, char **argv)
{
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  // define some Epetra objects

  int n = Comm.NumProc() * 4;
  Epetra_Map Map(n, 0, Comm);
  Epetra_MultiVector x(Map, 2); x.Random();
  Epetra_MultiVector b(Map, 2); x.Random();
  Epetra_CrsMatrix Matrix(Copy, Map, 0);
  // diagonal matrix
  for (int i = 0; i < Map.NumMyElements(); ++i)
  {
    int ii = Map.GID(i);
    double one = 1.0;
    Matrix.InsertGlobalValues(ii, 1, &one, &ii);
  }
  Matrix.FillComplete();

  Teuchos::ParameterList List;
  List.set("int parameter", 10);
  List.set("double parameter", 10.0);
  List.set("std::string parameter", "std::string");

  // ========================= //
  // Part I: generate XML file //
  // ========================= //
  
  EpetraExt::XMLWriter XMLWriter(Comm, "data.xml");

  std::vector<std::string> Content;
  Content.push_back("This is an example of description");
  Content.push_back("The description is as long as desired,");
  Content.push_back("just put it in a std::vector of strings.");

  XMLWriter.Create("MyProblem");
  XMLWriter.Write("Author", "myself and others");
  XMLWriter.Write("Date", "May 2006");
  XMLWriter.Write("MyMap", Map);
  XMLWriter.Write("MyMatrix", Matrix);
  XMLWriter.Write("MyLHS", x);
  XMLWriter.Write("MyRHS", b);
  XMLWriter.Write("MyContent", Content);
  XMLWriter.Write("MyParameters", List);
  XMLWriter.Close();

  // ================== //
  // Part II: read data //
  // ================== //
  
  EpetraExt::XMLReader XMLReader(Comm, "data.xml");

  Epetra_Map* MyMap;
  Epetra_CrsMatrix* MyMatrix;
  Epetra_MultiVector* MyLHS;
  Epetra_MultiVector* MyRHS;
  Teuchos::ParameterList MyParameters;
  std::vector<std::string> Author;
  std::vector<std::string> Date;
  std::vector<std::string> MyContent;

  XMLReader.Read("Author", Author);
  XMLReader.Read("Date", Date);
  XMLReader.Read("MyMap", MyMap);
  XMLReader.Read("MyMatrix", MyMatrix);
  XMLReader.Read("MyLHS", MyLHS);
  XMLReader.Read("MyRHS", MyRHS);
  XMLReader.Read("MyContent", MyContent);
  XMLReader.Read("MyParameters", MyParameters);

  std::cout << *MyMap;
  std::cout << *MyMatrix;
  std::cout << *MyLHS;
  std::cout << *MyRHS;
  if (Comm.MyPID() == 0)
  {
    int Msize = (int) MyContent.size();
    for (int i = 0; i < Msize; ++i)
      std::cout << MyContent[i] << std::endl;

    std::cout << MyParameters;
    std::cout << "Author = " << Author[0] << std::endl;
    std::cout << "Date   = " << Date[0] << std::endl;
  }

  delete MyMap;
  delete MyMatrix;
  delete MyLHS;
  delete MyRHS;

#ifdef HAVE_MPI
  MPI_Finalize();
#endif

//.........这里部分代码省略.........

void Filter::Kalman(Joint joint, double &dx, double &dy)
{
	Kalmans[Kalman_count++] = joint;
	Kalman_count = Kalman_count % Kalman_limit;
	Kalman_num++;
	if (Kalman_num > Kalman_limit)
	{
		Kalman_num = Kalman_limit;
	}
	if (Kalman_num < Kalman_limit)
	{
		dx = joint.Position.X;
		dy = joint.Position.Y;
		return;
	}
	else
	{
		//X, Y
		int haha;
		haha = 1;
		double x[5], y[5];
		int pos = Kalman_count;
		for (int i = 0; i < 5; i++)
		{
			x[i] = Kalmans[pos].Position.X;
			y[i] = Kalmans[pos].Position.Y;
			pos++;
			pos = pos%Kalman_limit;
		}
		//求系数Ax, Ay
		double Ax[5] = {
			/*a0*/ x[0],
			/*a1*/ 4 * (x[1] - x[0]) - 3 * x[2] + 4 * x[3] / 3 - x[4] / 4,
			/*a2*/ 11 * x[4] / 24 - 7 * x[3] / 3 + 19 * x[2] / 4 - 13 * (x[1] - x[0]) / 3,
			/*a3*/ x[4] / 3 - 7 * x[3] / 6 + x[2] - (x[1] - x[0]) / 2,
			/*a4*/ (x[4] - 4 * x[3] + 6 * x[2] + 4 * (x[1] - x[0])) / 24
		};
		double Ay[5] = {
			/*a0*/ y[0],
			/*a1*/ 4 * (y[1] - y[0]) - 3 * y[2] + 4 * y[3] / 3 - y[4] / 4,
			/*a2*/ 11 * y[4] / 24 - 7 * y[3] / 3 + 19 * y[2] / 4 - 13 * (y[1] - y[0]) / 3,
			/*a3*/ y[4] / 3 - 7 * y[3] / 6 + y[2] - (y[1] - y[0]) / 2,
			/*a4*/ (y[4] - 4 * y[3] + 6 * y[2] + 4 * (y[1] - y[0])) / 24
		};

		//求转换矩阵Fx, Fy
		Matrix Fx(4, 4,
							 new double[16]{
			  1, 1, -0.5, (Ax[1] + 6 * Ax[3] - 4 * Ax[4]) / (24 * Ax[4]),
				0, 1, 1, 0.5,
				0, 0, 1, 1,
				0, 0, 0, 1});
		Matrix Fy(4, 4,
							 new double[16]{
			1, 1, -0.5, (Ay[1] + 6 * Ay[3] - 4 * Ay[4]) / (24 * Ay[4]),
				0, 1, 1, 0.5,
				0, 0, 1, 1,
				0, 0, 0, 1});
		//求ε(t|t-1)
		Matrix ex(4, 4), ey(4, 4);
		ex = Fx*Kalman_ex*(!Fx);
		ey = Fy*Kalman_ey*(!Fy);
		//cout << "ex" << endl; ex.print();
		//cout << "ey" << endl; ey.print();

		Matrix Bx(4, 1), By(4, 1);
		//cout << "!Kalman_C" << endl; (!Kalman_C).print();
		//cout << "Kalman_vx" << endl; Kalman_vx.print();
		//cout << "Kalman_C" << endl; Kalman_C.print();
		//cout << "!Kalman_C" << endl; (!Kalman_C).print();
		//cout << "Kalman_C*ex" << endl; (Kalman_C*ex).print();
		//cout << "Kalman_C*ex*(!Kalman_C)" << endl; (Kalman_C*ex*(!Kalman_C)).print();
		//cout << "(~(Kalman_vx + Kalman_C*ex*(!Kalman_C)))" << endl;
		//(~(Kalman_vx + Kalman_C*ex*(!Kalman_C))).print();
		Bx = ex*(!Kalman_C)*(~(Kalman_vx + Kalman_C*ex*(!Kalman_C)));
		//cout << "Bx" << endl; Bx.print();
		By = ey*(!Kalman_C)*(~(Kalman_vy + Kalman_C*ey*(!Kalman_C)));
		
		Matrix I4(4, 4);
		I4.SetIdentity();
		Kalman_Sx = (I4 - Bx*Kalman_C)*(Fx*Kalman_Sx + Kalman_Gx) +
			Bx*Matrix(1, 1, new double[1] {joint.Position.X});
		//cout << "Kalman_Sx" << endl; Kalman_Sx.print();
		Kalman_Sy = (I4 - By*Kalman_C)*(Fy*Kalman_Sy + Kalman_Gy) +
			By*Matrix(1, 1, new double[1] {joint.Position.Y});

		Kalman_ex = ex - Bx*Kalman_C*ex;
		Kalman_ey = ey - By*Kalman_C*ey;

		dx = Kalman_Sx.at(0, 0);
		dy = Kalman_Sy.at(0, 0);
	}
}

inline void  GraphicsSystem::MatrixPush()
{
    MatrixStack << Matrix(MatrixStack[curMatrix]);
    ++curMatrix;
}

//==========================================================================
int Ifpack_ICT::Compute() 
{
  if (!IsInitialized()) 
    IFPACK_CHK_ERR(Initialize());

  Time_.ResetStartTime();
  IsComputed_ = false;

  NumMyRows_ = A_.NumMyRows();
  int Length = A_.MaxNumEntries();
  vector<int>    RowIndices(Length);
  vector<double> RowValues(Length);

  bool distributed = (Comm().NumProc() > 1)?true:false;

  if (distributed)
  {
    SerialComm_ = Teuchos::rcp(new Epetra_SerialComm);
    SerialMap_ = Teuchos::rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
    assert (SerialComm_.get() != 0);
    assert (SerialMap_.get() != 0);
  }
  else
    SerialMap_ = Teuchos::rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);

  int RowNnz;
#ifdef IFPACK_FLOPCOUNTERS
  double flops = 0.0;
#endif

  H_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy,*SerialMap_,0));
  if (H_.get() == 0)
    IFPACK_CHK_ERR(-5); // memory allocation error

  // get A(0,0) element and insert it (after sqrt)
  IFPACK_CHK_ERR(A_.ExtractMyRowCopy(0,Length,RowNnz,
                                     &RowValues[0],&RowIndices[0]));

  // skip off-processor elements
  if (distributed)
  {
    int count = 0;
    for (int i = 0 ;i < RowNnz ; ++i) 
    {
      if (RowIndices[i] < NumMyRows_){
        RowIndices[count] = RowIndices[i];
        RowValues[count] = RowValues[i];
        ++count;
      }
      else
        continue;
    }
    RowNnz = count;
  }

  // modify diagonal
  double diag_val = 0.0;
  for (int i = 0 ;i < RowNnz ; ++i) {
    if (RowIndices[i] == 0) {
      double& v = RowValues[i];
      diag_val = AbsoluteThreshold() * EPETRA_SGN(v) +
        RelativeThreshold() * v;
      break;
    }
  }

  diag_val = sqrt(diag_val);
  int diag_idx = 0;
  EPETRA_CHK_ERR(H_->InsertGlobalValues(0,1,&diag_val, &diag_idx));

  int oldSize = RowNnz;

  // The 10 is just a small constant to limit collisons as the actual keys
  // we store are the indices and not integers
  // [0..A_.MaxNumEntries()*LevelofFill()].
  Ifpack_HashTable Hash( 10 * A_.MaxNumEntries() * LevelOfFill(), 1);

  // start factorization for line 1
  for (int row_i = 1 ; row_i < NumMyRows_ ; ++row_i) {

    // get row `row_i' of the matrix
    IFPACK_CHK_ERR(A_.ExtractMyRowCopy(row_i,Length,RowNnz,
                                       &RowValues[0],&RowIndices[0]));

    // skip off-processor elements
    if (distributed)
    {
      int count = 0;
      for (int i = 0 ;i < RowNnz ; ++i) 
      {
        if (RowIndices[i] < NumMyRows_){
          RowIndices[count] = RowIndices[i];
          RowValues[count] = RowValues[i];
          ++count;
        }
        else
          continue;
      }
      RowNnz = count;
//.........这里部分代码省略.........

//==============================================================================
int Ifpack_PointRelaxation::
ApplyInverseSGS_RowMatrix(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{
  int NumVectors = X.NumVectors();
  int Length = Matrix().MaxNumEntries();
  vector<int> Indices(Length);
  vector<double> Values(Length);

  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
  if (IsParallel_) {
    Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
  }
  else
    Y2 = Teuchos::rcp( &Y, false );

  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
  X.ExtractView(&x_ptr);
  Y.ExtractView(&y_ptr);
  Y2->ExtractView(&y2_ptr);
  Diagonal_->ExtractView(&d_ptr);
  
  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
    
    // only one data exchange per sweep
    if (IsParallel_)
      IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));

    for (int i = 0 ; i < NumMyRows_ ; ++i) {

      int NumEntries;
      int col;
      double diag = d_ptr[i];

      IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
                                               &Values[0], &Indices[0]));

      for (int m = 0 ; m < NumVectors ; ++m) {

        double dtemp = 0.0;

        for (int k = 0 ; k < NumEntries ; ++k) {

          col = Indices[k];
          dtemp += Values[k] * y2_ptr[m][col];
        }

        y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
      }
    }

    for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {

      int NumEntries;
      int col;
      double diag = d_ptr[i];

      IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
                                               &Values[0], &Indices[0]));

      for (int m = 0 ; m < NumVectors ; ++m) {

        double dtemp = 0.0;
        for (int k = 0 ; k < NumEntries ; ++k) {

          col = Indices[k];
          dtemp += Values[k] * y2_ptr[m][col];
        }

        y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;

      }
    }

    if (IsParallel_)
      for (int m = 0 ; m < NumVectors ; ++m) 
        for (int i = 0 ; i < NumMyRows_ ; ++i)
          y_ptr[m][i] = y2_ptr[m][i];
  }

  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
  return(0);
}

returnValue CondensingExport::setupEvaluation( )
{
	x.setup( "x",        (getN()+1), getNX(), REAL,ACADO_VARIABLES );
	u.setup( "u",        getN(), getNU(),     REAL,ACADO_VARIABLES );
	p.setup( "p",        1, getNP(),          REAL,ACADO_VARIABLES );

	state.setup   ( "state",    1,getNX()*(getNX()+getNU()+1) + getNU() +getNP(), REAL,ACADO_WORKSPACE );

	if ( isInitialStateFixed( ) == BT_FALSE )
	{
		x0Ref.setup ( "x0Ref",  1, getNX(), REAL,ACADO_VARIABLES );
		x0Ref2.setup( "x0Ref2", 1, getNX(), REAL,ACADO_VARIABLES );
	}

	E.setup  ( "E",   getN()*getNX(), getN()*getNU(), REAL,ACADO_WORKSPACE );
	lbA.setup( "lbA", getNumStateBounds(), 1, REAL,ACADO_WORKSPACE );
	ubA.setup( "ubA", getNumStateBounds(), 1, REAL,ACADO_WORKSPACE );

	Matrix zeroXU = zeros( getNX(),getNU() );
	Matrix idX    = eye( getNX() );

	//
	// setupQP
	//
	setupQP.setup( "setupQP" );
    
	if ( performsSingleShooting() == BT_TRUE )
		setupQP.addStatement( state.getCols( 0,getNX() ) == x.getRow(0) );

	// TODO: this part should be preinitialized. More specifically, E & QE matrices should be preinitializes to 0.
	uint run1, run2;
	for( run1 = 0; run1 < getN()-1; run1++ )
		for( run2 = 1+run1; run2 < getN(); run2++ )
			setupQP.addStatement( E.getSubMatrix( run1*getNX(),(run1+1)*getNX(), run2*getNU(),(run2+1)*getNU() ) == zeroXU );


    // Write state bounds to the file
	// TODO: Since the bounds are fixed in this case, they should be preinitialized
	if( getNumStateBounds( ) > 0 )
	{
		Vector xLowerBounds(nxBounds), xUpperBounds(nxBounds);
		for( run1 = 0; run1 < nxBounds; run1++ )
		{
			xLowerBounds(run1) = xBounds.getLowerBound( xBoundsIdx[run1]/getNX(),xBoundsIdx[run1]%getNX() );
			xUpperBounds(run1) = xBounds.getUpperBound( xBoundsIdx[run1]/getNX(),xBoundsIdx[run1]%getNX() );
		}
		
		setupQP.addStatement( lbA == xLowerBounds );
		setupQP.addStatement( ubA == xUpperBounds );
	}
	setupQP.addLinebreak( );


	if ( isInitialStateFixed( ) == BT_FALSE )
	{
		setupQP.addStatement( Dx0  == x.getRow(0) - x0Ref );
		setupQP.addStatement( Dx0b == x.getRow(0) - x0Ref2 );
		setupQP.addLinebreak( );
	}
	
	
	// compute QQF if necessary
	if ( QQF.isGiven( ) == BT_FALSE )
		setupQP.addStatement( QQF == Q + QF );

	ExportIndex reset( String("reset") );
	setupQP.addIndex( reset );
	setupQP.addStatement( reset == 1 );

	ExportIndex run( String("run1") );
	ExportForLoop loop( run, 0,getN() );
	
	setupQP.addIndex( run );

	if ( performsSingleShooting() == BT_FALSE ) {
		loop.addStatement( reset == 1 );
		loop.addStatement( state.getCols( 0,getNX() ) == x.getRow( run ) );
	}
	
	// no free parameters implemented yet!
	uint uIdx = getNX() * ( 1+getNX() );
	uint pIdx = getNX() * ( 1+getNX()+getNU() );
	uIdx = pIdx;
	pIdx = pIdx + getNU();
	loop.addStatement( state.getCols( uIdx,pIdx ) == u.getRow( run ) );
	loop.addStatement( state.getCols( pIdx,pIdx+getNP() ) == p );
	loop.addLinebreak( );
	
	if ( integrator->equidistantControlGrid() )
	{
		loop.addFunctionCall( "integrate", state, reset.makeArgument() );
	}
	else
	{
		loop.addFunctionCall( "integrate", state, reset.makeArgument(), run.makeArgument() );
	}
	if( performsSingleShooting() ) loop.addStatement( reset == 0 );
	
	if ( performsSingleShooting() == BT_TRUE )
	{
//.........这里部分代码省略.........

//==============================================================================
int Ifpack_PointRelaxation::Compute()
{
  int ierr = 0;
  if (!IsInitialized())
    IFPACK_CHK_ERR(Initialize());

  Time_->ResetStartTime();

  // reset values
  IsComputed_ = false;
  Condest_ = -1.0;

  if (NumSweeps_ == 0) ierr = 1; // Warning: no sweeps performed.
  if (NumSweeps_ < 0)
    IFPACK_CHK_ERR(-2); // at least one application
  
  Diagonal_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );

  if (Diagonal_ == Teuchos::null)
    IFPACK_CHK_ERR(-5);

  IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*Diagonal_));

  // check diagonal elements, store the inverses, and verify that
  // no zeros are around. If an element is zero, then by default
  // its inverse is zero as well (that is, the row is ignored).
  for (int i = 0 ; i < NumMyRows_ ; ++i) {
    double& diag = (*Diagonal_)[i];
    if (IFPACK_ABS(diag) < MinDiagonalValue_)
      diag = MinDiagonalValue_;
    if (diag != 0.0)
      diag = 1.0 / diag;
  }
  ComputeFlops_ += NumMyRows_;

#if 0
  // some methods require the inverse of the diagonal, compute it
  // now and store it in `Diagonal_'
  if ((PrecType_ == IFPACK_JACOBI) || (PrecType_ == IFPACK_GS)) {
    Diagonal_->Reciprocal(*Diagonal_);
    // update flops
    ComputeFlops_ += NumMyRows_;
  }
#endif

  // We need to import data from external processors. Here I create an
  // Epetra_Import object because I cannot assume that Matrix_ has one.
  // This is a bit of waste of resources (but the code is more robust).
  // Note that I am doing some strange stuff to set the components of Y
  // from Y2 (to save some time).
  //
  if (IsParallel_ && ((PrecType_ == IFPACK_GS) || (PrecType_ == IFPACK_SGS))) {
    Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
                                  Matrix().RowMatrixRowMap()) );
    if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
  }

  ++NumCompute_;
  ComputeTime_ += Time_->ElapsedTime();
  IsComputed_ = true;

  IFPACK_CHK_ERR(ierr);
  return(0);
}

//==============================================================================
int Ifpack_PointRelaxation::
ApplyInverseGS_RowMatrix(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{
  int NumVectors = X.NumVectors();

  int Length = Matrix().MaxNumEntries();
  vector<int> Indices(Length);
  vector<double> Values(Length);

  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
  if (IsParallel_)
    Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
  else
    Y2 = Teuchos::rcp( &Y, false );

  // extract views (for nicer and faster code)
  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
  X.ExtractView(&x_ptr);
  Y.ExtractView(&y_ptr);
  Y2->ExtractView(&y2_ptr);
  Diagonal_->ExtractView(&d_ptr);

  for (int j = 0; j < NumSweeps_ ; j++) {

    // data exchange is here, once per sweep
    if (IsParallel_)
      IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));

    // FIXME: do I really need this code below?
    if (NumVectors == 1) {

      double* y0_ptr = y_ptr[0];
      double* y20_ptr = y2_ptr[0];
      double* x0_ptr = x_ptr[0];

      if(!DoBackwardGS_){      
        /* Forward Mode */
        for (int i = 0 ; i < NumMyRows_ ; ++i) {

          int NumEntries;
          int col;
          IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
                                                   &Values[0], &Indices[0]));
          
          double dtemp = 0.0;
          for (int k = 0 ; k < NumEntries ; ++k) {
            
            col = Indices[k];
            dtemp += Values[k] * y20_ptr[col];
          }
          
          y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
        }
      }
      else {
        /* Backward Mode */
        for (int i = NumMyRows_  - 1 ; i > -1 ; --i) {

          int NumEntries;
          int col;
          IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
                                                   &Values[0], &Indices[0]));
          double dtemp = 0.0;
          for (int k = 0 ; k < NumEntries ; ++k) {

            col = Indices[k];
            dtemp += Values[k] * y20_ptr[i];
          }
          
          y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
        }
      }
      
      // using Export() sounded quite expensive
      if (IsParallel_)
        for (int i = 0 ; i < NumMyRows_ ; ++i)
          y0_ptr[i] = y20_ptr[i];

    }
    else {
      if(!DoBackwardGS_){      
        /* Forward Mode */
        for (int i = 0 ; i < NumMyRows_ ; ++i) {
          
          int NumEntries;
          int col;
          IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
                                                   &Values[0], &Indices[0]));
          
          for (int m = 0 ; m < NumVectors ; ++m) {
            
            double dtemp = 0.0;
            for (int k = 0 ; k < NumEntries ; ++k) {
              
              col = Indices[k];
              dtemp += Values[k] * y2_ptr[m][col];
            }
            
            y2_ptr[m][i] += DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
//.........这里部分代码省略.........

//==========================================================================
int Ifpack_SPARSKIT::Compute() 
{
  if (!IsInitialized()) 
    IFPACK_CHK_ERR(Initialize());

  IsComputed_ = false;

  // convert the matrix into SPARSKIT format. The matrix is then
  // free'd after method Compute() returns.

  // convert the matrix into CSR format. Note that nnz is an over-estimate,
  // since it contains the nonzeros corresponding to external nodes as well.
  int n   = Matrix().NumMyRows();
  int nnz = Matrix().NumMyNonzeros();

  vector<double> a(nnz);
  vector<int>    ja(nnz);
  vector<int>    ia(n + 1);

  const int MaxNumEntries = Matrix().MaxNumEntries();

  vector<double> Values(MaxNumEntries);
  vector<int>    Indices(MaxNumEntries);

  int count = 0;

  ia[0] = 1;
  for (int i = 0 ; i < n ; ++i)
  {
    int NumEntries;
    int NumMyEntries = 0;
    Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0], 
                              &Indices[0]);

    // NOTE: There might be some issues here with the ILU(0) if the column indices aren't sorted.
    // The other factorizations are probably OK.

    for (int j = 0 ; j < NumEntries ; ++j)
    {
      if (Indices[j] < n) // skip non-local columns
      {
        a[count]  = Values[j];
        ja[count] = Indices[j] + 1; // SPARSKIT is FORTRAN
        ++count;
        ++NumMyEntries;
      }
    }
    ia[i + 1] = ia[i] + NumMyEntries;
  }

  if (mbloc_ == -1) mbloc_ = n;

  int iwk;

  if (Type_ == "ILUT" || Type_ == "ILUTP" || Type_ == "ILUD" ||
      Type_ == "ILUDP")
    iwk = nnz + 2 * lfil_ * n;
  else if (Type_ == "ILUK")
    iwk = (2 * lfil_ + 1) * nnz + 1;
  else if (Type_ == "ILU0")
    iwk = nnz+2;

  int ierr = 0;

  alu_.resize(iwk);
  jlu_.resize(iwk);
  ju_.resize(n + 1);

  vector<int>    jw(n + 1);
  vector<double> w(n + 1);

  if (Type_ == "ILUT")
  {
    jw.resize(2 * n);
    F77_ILUT(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_,
             &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
  }
  else if (Type_ == "ILUD")
  {
    jw.resize(2 * n);
    F77_ILUD(&n, &a[0], &ja[0], &ia[0], &alph_, &tol_,
             &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
  }
  else if (Type_ == "ILUTP")
  {
    jw.resize(2 * n);
    iperm_.resize(2 * n);
    F77_ILUTP(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_, &permtol_, 
              &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0],
              &iperm_[0], &ierr);
    for (int i = 0 ; i < n ; ++i)
      iperm_[i]--;
  }
  else if (Type_ == "ILUDP")
  {
    jw.resize(2 * n);
    iperm_.resize(2 * n);
    F77_ILUDP(&n, &a[0], &ja[0], &ia[0], &alph_, &droptol_, &permtol_, 
              &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &n, &w[0], &jw[0],
//.........这里部分代码省略.........

    FwdPeriodAdapter::FwdPeriodAdapter(
                               const ext::shared_ptr<MarketModel>& largeModel,
                               Size period,
                               Size offset,
                               const std::vector<Spread>& newDisplacements)
    :
      numberOfFactors_(largeModel->numberOfFactors()),
          numberOfRates_((largeModel->numberOfRates()-offset) / (period > 0 ? period : 1) ),
      numberOfSteps_(largeModel->numberOfSteps()),
      pseudoRoots_(numberOfSteps_, Matrix(numberOfRates_,
                                          numberOfFactors_)),
                                          displacements_(newDisplacements)
    {
        QL_REQUIRE( period >0, "period must  be greater than zero in fwdperiodadapter");
        QL_REQUIRE(period > offset, "period must be greater than offset in fwdperiodadapter");

        const std::vector<Spread>& largeDisplacements_ =
            largeModel->displacements();

        if (displacements_.size() == 1)
        {
            Real dis = displacements_[0];
            displacements_.resize(numberOfRates_);
            std::fill(displacements_.begin(), displacements_.end(), dis);
        }

        if (displacements_.size() ==0) // if not specified use average across rate
        {
            displacements_.reserve(numberOfRates_);
            Size m=0;
            Real sum=0.0;
            for (Size k=0; k < numberOfRates_; ++k)
            {
                for (Size l=0; l < period; ++l, ++m)
                    sum+= largeDisplacements_[m];

                displacements_.push_back(sum/period);
            }
        }
        QL_REQUIRE( displacements_.size() == numberOfRates_,"newDisplacements should be empty,1, or number of new rates in fwdperiodadapter");

        LMMCurveState largeCS(largeModel->evolution().rateTimes());
        largeCS.setOnForwardRates(largeModel->initialRates());

        LMMCurveState smallCS(
                ForwardForwardMappings::RestrictCurveState(largeCS,
                                    period, offset
                                        ));

        initialRates_ =smallCS.forwardRates();

        Real finalReset = smallCS.rateTimes()[smallCS.numberOfRates()-1];
        std::vector<Time> oldEvolutionTimes(largeModel->evolution().evolutionTimes());
        std::vector<Time> newEvolutionTimes;
        for (Size i =0; i < oldEvolutionTimes.size() && oldEvolutionTimes[i]<= finalReset; ++i)
            newEvolutionTimes.push_back(oldEvolutionTimes[i]);

        evolution_=EvolutionDescription(smallCS.rateTimes(),
                                        newEvolutionTimes);

        numberOfSteps_ = newEvolutionTimes.size();


        const std::vector<Time>& rateTimes =
            smallCS.rateTimes();
        // we must ensure we step through all rateTimes
        const std::vector<Time>& evolutionTimes =
            evolution_.evolutionTimes();

        std::set<Time> setTimes(evolutionTimes.begin(),evolutionTimes.end());

        for (Size i=0; i < rateTimes.size()-1; ++i)
            QL_REQUIRE(setTimes.find(rateTimes[i]) != setTimes.end(),
                        "every new rate time except last must be an evolution time in fwdperiod adapter");


        Matrix YMatrix =
            ForwardForwardMappings::YMatrix( largeCS,
                                                 largeDisplacements_,
                                                  displacements_,
                                                 period,
                                                 offset
                                                 );

        const std::vector<Size>& alive =
            evolution_.firstAliveRate();

        for (Size k = 0; k<numberOfSteps_; ++k) {
            pseudoRoots_[k]=YMatrix*largeModel->pseudoRoot(k);
            for (Size i=0; i<alive[k]; ++i)
                std::fill(pseudoRoots_[k].row_begin(i),
                          pseudoRoots_[k].row_end(i),
                          0.0);
        }
    }

std::vector<double> NNTrainer::gradient( const double lambda)
{
    int numExamples = trainingSet->size();
    //-- Create n matrices with the dimensions of the weight matrices:
    std::vector<Matrix> Delta;

    for (int i = 0; i < (int) nn->getWeights().size(); i++)
	Delta.push_back( Matrix( nn->getWeights().at(i)->getNumRows(), nn->getWeights().at(i)->getNumCols() ));

    //-- Iterate over all training examples
    for (int i = 0; i < numExamples; i++)
    {
	//-- Forward-propagate the network:
	nn->setInput( trainingSet->at(i).x );

	//-- Create vector to store the increments:
	//-- Increments will be stored in reverse order (i.e. the last increment first)
	std::vector<Matrix> inc;

	//-- Increment for output layer
	Matrix output = Matrix(nn->getOutput(), nn->getOutput().size(), 1);
	Matrix y = Matrix(trainingSet->at(i).y , trainingSet->at(i).y.size(), 1);

	inc.push_back( output - y);

	//-- Increment for hidden layers
	for (int l = nn->getL() - 2; l > 0; l--)
	{

	    Matrix aux1 = nn->getWeights().at(l)->transpose() * inc.back();
	    Matrix aux2( aux1.getNumRows()-1, aux1.getNumCols());

	    for (int j = 0; j < aux2.getNumCols(); j++)
		for (int k = 0; k < aux2.getNumRows(); k++)
		    aux2.set( k, j, aux1.get(k+1, j) * sigmoidGradient( nn->getActivation(l).at(k)) );

	    inc.push_back( aux2 );
	}

	//-- Input layer has no error associated (has no weights associated)

	//-- Accumulate error:
	for (int l = 0; l < (int) Delta.size(); l++)
	{
	    Matrix aux1( Delta.at(l).getNumRows(), Delta.at(l).getNumCols() );

	    for (int j = 0; j < aux1.getNumRows(); j++)
		aux1.set( j, 0, inc.at( inc.size()- l -1).get( j, 0) );

	    for (int j = 0; j < aux1.getNumRows(); j++)
		for (int k = 1; k < aux1.getNumCols(); k++)
		    aux1.set( j, k, inc.at( inc.size()- l -1).get( j, 0) * nn->getActivation(l).at(k-1));

	    Delta.at(l) += aux1;
	}


    }

    //-- Divide by number of training examples:
    for (int l = 0; l < (int) Delta.size(); l++)
	Delta.at(l) /= numExamples;

    //-- Regularization
    //------------------------------------------------------------------------
    if (lambda != 0)
    {
	for (int l = 0; l < (int) Delta.size(); l++)
	{
	    Matrix aux(nn->getWeights().at(l)->getNumRows(), nn->getWeights().at(l)->getNumCols() );

	    for (int j = 0; j < aux.getNumRows(); j++)
		for (int k = 1; k < aux.getNumCols(); k++)
		    aux.set( j, k, nn->getWeights().at(l)->get(j, k) * lambda / numExamples);

	    Delta.at(l) += aux;
	}
    }


    //-- Unroll gradient:
    //---------------------------------------------------------------------------
    std::vector<double> unrolled = Delta.front().unroll();

    for (int l = 1; l < (int) Delta.size(); l++)
	for (int j = 0; j < Delta.at(l).getNumRows(); j++)
	    for (int k = 0; k < Delta.at(l).getNumCols(); k++)
		unrolled.push_back( Delta.at(l).get(j, k));


    return unrolled;
}