[Trilinos-Users] Epetra_RowMatrix

[Mitchell, John A]

Hello,

I am having trouble with a "scope" issue as it relates to creating an Epetra_RowMatrix.   

There is an enclosing "block" (brackets) that constructs an Epetra_FEVbrMatrix*

Inside the "block" I fill the matrix with numerical values and then extract and print the diagonal.  OK so far.

Immediately after exiting the block, I extract and print the diagonal again but for some reason it is now zero?  

I'm missing something here?  

Thanks for your help.
John



	Epetra_RowMatrix *mPtr;

	{
		const Epetra_BlockMap& rowMap   = jacobian->getRowMap();
		const Epetra_BlockMap& colMap   = jacobian->getColMap();

		/*
		 * Epetra Matrix
		 * PERHAPS the 'operator' can keep its own copy of Graph, then this
		 * can be a 'View'
		 */

		Epetra_FEVbrMatrix *m = new Epetra_FEVbrMatrix(Copy,*(graphPtr.get()));

		Epetra_SerialDenseMatrix k;
		k.Shape(vectorNDF,vectorNDF);
		for(int row=0;row<rowMap.NumMyElements();row++){
			Pd_shared_ptr_Array<int> rowLIDs = jacobian->getColumnLIDs(row);
			std::size_t numCol = rowLIDs.getSize();
			int *cols = rowLIDs.get();
			BOOST_CHECK(0==m->BeginReplaceMyValues(row,numCol,cols));
			/*
			 * loop over columns in row and submit block entry
			 */
			const Pd_shared_ptr_Array<double>& actualK = jacobian->computeRowStiffness(row, rowLIDs);
			const double *kPtr = actualK.get();

			for(std::size_t col=0;col<numCol;col++){
				/*
				 * Fill 'k'
				 */
				for(int c=0;c<3;c++){
					double *colK = k[c];
					for(int r=0;r<3;r++,kPtr++)
						colK[r] = *kPtr;
				}
				BOOST_CHECK(0==m->SubmitBlockEntry(k));
			}
			/*
			 * Finalize this row
			 */
			BOOST_CHECK(0==m->EndSubmitEntries());
		}

		BOOST_CHECK(0==m->FillComplete());
		mPtr = m;

		Epetra_Vector stiffnessDiag(colMap,true);
		mPtr->ExtractDiagonalCopy(stiffnessDiag);
		BOOST_CHECK(stiffnessDiag.GlobalLength()==numPoints*3);
		std::cout << "A -- diag\n";
		std::cout << "Stiffness Diagonal \n";
		for(int n=0;n<stiffnessDiag.GlobalLength();n++){
			std::cout << "n = " << n << "; diag[n] = " << stiffnessDiag[n] << "\n";
		}


	}

	/*
	 * extract and print diagonal again
	 */
	const Epetra_Map& colMap = mPtr->RowMatrixColMap();
	Epetra_Vector stiffnessDiag(colMap,true);
	mPtr->ExtractDiagonalCopy(stiffnessDiag);
	BOOST_CHECK(stiffnessDiag.GlobalLength()==numPoints*3);
	std::cout << "B -- diag\n";
	std::cout << "Stiffness Diagonal \n";
	for(int n=0;n<stiffnessDiag.GlobalLength();n++){
		std::cout << "n = " << n << "; diag[n] = " << stiffnessDiag[n] << "\n";
	}