C++ DataVector::GetRows方法代码示例

void SolveEvolutionMatrix(
	DataMatrix<double> & matM,
	DataMatrix<double> & matB,
	DataVector<double> & vecAlphaR,
	DataVector<double> & vecAlphaI,
	DataVector<double> & vecBeta,
	DataMatrix<double> & matVR
) {
	int iInfo = 0;

	char jobvl = 'N';
	char jobvr = 'V';

	int n = matM.GetRows();
	int lda = matM.GetRows();
	int ldb = matM.GetRows();
	int ldvl = 1;
	int ldvr = matVR.GetRows();

	DataVector<double> vecWork;
	vecWork.Initialize(8 * n);
	int lwork = vecWork.GetRows();

	dggev_(
		&jobvl,
		&jobvr,
		&n,
		&(matM[0][0]),
		&lda,
		&(matB[0][0]),
		&ldb,
		&(vecAlphaR[0]),
		&(vecAlphaI[0]),
		&(vecBeta[0]),
		NULL,
		&ldvl,
		&(matVR[0][0]),
		&ldvr,
		&(vecWork[0]),
		&lwork,
		&iInfo);

	int nCount = 0;
	for (int i = 0; i < n; i++) {
		if (vecBeta[i] != 0.0) {
			//printf("%i %1.5e %1.5e\n", i, vecAlphaR[i] / vecBeta[i], vecAlphaI[i] / vecBeta[i]);
			nCount++;
		}
	}
/*
	for (int i = 0; i < 40; i++) {
		printf("%1.5e %1.5e\n", matVR[11][4*i+2], matVR[12][4*i+2]);
	}
*/
	Announce("%i total eigenvalues found", nCount);
}

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

	DataVector<double> dWI;
	dWI.Initialize(param.nPhiElements);

	DataVector<double> dRhoR;
	dRhoR.Initialize(param.nPhiElements);

	DataVector<double> dRhoI;
	dRhoI.Initialize(param.nPhiElements);

	// Loop over all horizontal wave numbers
	for (int nK = nKmin; nK <= nKmax; nK++) {

		// Build matrices
		char szMessage[100];
		sprintf(szMessage, "Building evolution matrices (k = %i)", nK);
		AnnounceStartBlock(szMessage);

		DataMatrix<double> matM;
		DataMatrix<double> matB;

		GenerateEvolutionMatrix(nK, param, matM, matB);

		AnnounceEndBlock("Done");

		// Solve the matrices
		AnnounceStartBlock("Solving evolution matrices");

		DataVector<double> vecAlphaR;
		DataVector<double> vecAlphaI;
		DataVector<double> vecBeta;
		DataMatrix<double> matVR;

		vecAlphaR.Initialize(matM.GetRows());
		vecAlphaI.Initialize(matM.GetRows());
		vecBeta  .Initialize(matM.GetRows());
		matVR    .Initialize(matM.GetRows(), matM.GetColumns());

		SolveEvolutionMatrix(
			matM,
			matB,
			vecAlphaR,
			vecAlphaI,
			vecBeta,
			matVR);

		// Sort eigenvalues
		std::multimap<double, int> mapSortedRows;
		for (int i = 0; i < vecBeta.GetRows(); i++) {
			if (vecBeta[i] != 0.0) {

				double dLambdaR = vecAlphaR[i] / vecBeta[i];
				double dLambdaI = vecAlphaI[i] / vecBeta[i];

				double dMR = dLambdaI;
				double dMI = - dLambdaR - 1.0;

				double dEvalMagnitude = fabs(dMR);

				//printf("\n%1.5e %1.5e ", dMR, dMI);

				// Purely imaginary eigenvalue
				if (dMR == 0.0) {
					// Wave must decay with height
					if (dMI < 0.0) {
						continue;

void GridPatchCartesianGLL::InterpolateData(
	const DataVector<double> & dAlpha,
	const DataVector<double> & dBeta,
	const DataVector<int> & iPatch,
	DataType eDataType,
	DataLocation eDataLocation,
	bool fInterpAllVariables,
	DataMatrix3D<double> & dInterpData,
	bool fIncludeReferenceState,
	bool fConvertToPrimitive
) {
	if (dAlpha.GetRows() != dBeta.GetRows()) {
		_EXCEPTIONT("Point vectors must have equivalent length.");
	}

	// Vector for storage interpolated points
	DataVector<double> dAInterpCoeffs;
	dAInterpCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dBInterpCoeffs;
	dBInterpCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dADiffCoeffs;
	dADiffCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dBDiffCoeffs;
	dBDiffCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dAInterpPt;
	dAInterpPt.Initialize(m_nHorizontalOrder);

	// Element-wise grid spacing
	double dDeltaA = GetElementDeltaA();
	double dDeltaB = GetElementDeltaB();

	// Physical constants
	const PhysicalConstants & phys = m_grid.GetModel().GetPhysicalConstants();

	// Loop throught all points
	for (int i = 0; i < dAlpha.GetRows(); i++) {

		// Element index
		if (iPatch[i] != GetPatchIndex()) {
			continue;
		}

		// Verify point lies within domain of patch
		const double Eps = 1.0e-10;
		if ((dAlpha[i] < m_box.GetAEdge(m_box.GetAInteriorBegin()) - Eps) ||
			(dAlpha[i] > m_box.GetAEdge(m_box.GetAInteriorEnd()) + Eps) ||
			(dBeta[i] < m_box.GetBEdge(m_box.GetBInteriorBegin()) - Eps) ||
			(dBeta[i] > m_box.GetBEdge(m_box.GetBInteriorEnd()) + Eps)
		) {
			_EXCEPTIONT("Point out of range");
		}

		// Determine finite element index
		int iA =
			(dAlpha[i] - m_box.GetAEdge(m_box.GetAInteriorBegin())) / dDeltaA;

		int iB =
			(dBeta[i] - m_box.GetBEdge(m_box.GetBInteriorBegin())) / dDeltaB;

		// Bound the index within the element
		if (iA < 0) {
			iA = 0;
		}
		if (iA >= (m_box.GetAInteriorWidth() / m_nHorizontalOrder)) {
			iA = m_box.GetAInteriorWidth() / m_nHorizontalOrder - 1;
		}
		if (iB < 0) {
			iB = 0;
		}
		if (iB >= (m_box.GetBInteriorWidth() / m_nHorizontalOrder)) {
			iB = m_box.GetBInteriorWidth() / m_nHorizontalOrder - 1;
		}

		iA = m_box.GetHaloElements() + iA * m_nHorizontalOrder;
		iB = m_box.GetHaloElements() + iB * m_nHorizontalOrder;

		// Compute interpolation coefficients
		PolynomialInterp::LagrangianPolynomialCoeffs(
			m_nHorizontalOrder,
			&(m_box.GetAEdges()[iA]),
			dAInterpCoeffs,
			dAlpha[i]);

		PolynomialInterp::LagrangianPolynomialCoeffs(
			m_nHorizontalOrder,
			&(m_box.GetBEdges()[iB]),
			dBInterpCoeffs,
			dBeta[i]);

		int nComponents;
		int nRElements = m_grid.GetRElements();

		double ** pData2D;

		// State Data: Perform interpolation on all variables
		if (eDataType == DataType_State) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	const EquationSet & eqns = m_grid.GetModel().GetEquationSet();

	int nComponents = eqns.GetComponents();
	int nTracers = eqns.GetTracers();

	DataVector<double> dPointwiseState;
	dPointwiseState.Initialize(nComponents);

	DataVector<double> dPointwiseRefState;
	dPointwiseRefState.Initialize(nComponents);

	DataVector<double> dPointwiseTracers;
	if (m_datavecTracers.size() > 0) {
		dPointwiseTracers.Initialize(nTracers);
	}

	// Evaluate the state on model levels
	for (int k = 0; k < m_grid.GetRElements(); k++) {
	for (int i = 0; i < m_box.GetATotalWidth(); i++) {
	for (int j = 0; j < m_box.GetBTotalWidth(); j++) {

		// Evaluate pointwise state
		test.EvaluatePointwiseState(
			m_grid.GetModel().GetPhysicalConstants(),
			time,
			m_dataZLevels[k][i][j],
			m_dataLon[i][j],
			m_dataLat[i][j],
			dPointwiseState,
			dPointwiseTracers);

		eqns.ConvertComponents(phys, dPointwiseState);

		for (int c = 0; c < dPointwiseState.GetRows(); c++) {
			m_datavecStateNode[iDataIndex][c][k][i][j] = dPointwiseState[c];
		}

		// Evaluate reference state
		if (m_grid.HasReferenceState()) {
			test.EvaluateReferenceState(
				m_grid.GetModel().GetPhysicalConstants(),
				m_dataZLevels[k][i][j],
				m_dataLon[i][j],
				m_dataLat[i][j],
				dPointwiseRefState);

			eqns.ConvertComponents(phys, dPointwiseRefState);

			for (int c = 0; c < dPointwiseState.GetRows(); c++) {
				m_dataRefStateNode[c][k][i][j] = dPointwiseRefState[c];
			}
		}

		// Evaluate tracers
		for (int c = 0; c < dPointwiseTracers.GetRows(); c++) {
			m_datavecTracers[iDataIndex][c][k][i][j] = dPointwiseTracers[c];
		}
	}
	}
	}

	// Evaluate the state on model interfaces
	for (int k = 0; k <= m_grid.GetRElements(); k++) {
	for (int i = 0; i < m_box.GetATotalWidth(); i++) {
	for (int j = 0; j < m_box.GetBTotalWidth(); j++) {