[Trilinos-Users] Epetra_RowMatrix

Mitchell, John A jamitch at sandia.gov
Fri Apr 30 15:33:43 MDT 2010


Thanks Chris -- that was the issue.  

________________________________________
From: Baker, Christopher G. [bakercg at ornl.gov]
Sent: Friday, April 30, 2010 2:40 PM
To: Mitchell, John A; trilinos-users at software.sandia.gov
Subject: Re: [Trilinos-Users] Epetra_RowMatrix

Hi John,

The included code does not check the return value of ExtractDiagonalCopy() in either case, but more significantly, in the latter case. This means that an error in the routine will go unreported. In particular, I am concerned that the error is that the row map for the second/outer Epetra_Vector stiffnessDiag does not match that of the matrix; (in fact, the stiffnessDiag is being created from the column map, which is incorrect if the column map and the row map of the matrix differ.) If this is indeed the problem, then ExtractDiagonalCopy() will return an error -2.

Chris


On 4/30/10 4:29 PM, "Mitchell, John A" <jamitch at sandia.gov> wrote:

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

_______________________________________________
Trilinos-Users mailing list
Trilinos-Users at software.sandia.gov
http://software.sandia.gov/mailman/listinfo/trilinos-users






More information about the Trilinos-Users mailing list