FEI Version of the Day
test_Matrix.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 #include <fei_macros.hpp>
00044 #include <cmath>
00045 #include <fei_mpi.h>
00046 #include <test_utils/fei_test_utils.hpp>
00047 #include <test_utils/LibraryFactory.hpp>
00048 #include <test_utils/test_Matrix.hpp>
00049 #include <test_utils/test_VectorSpace.hpp>
00050 #include <test_utils/test_MatrixGraph.hpp>
00051 #include <fei_MatrixGraph_Impl2.hpp>
00052 #include <fei_Factory.hpp>
00053 #include <fei_defs.h>
00054 #include <snl_fei_Factory.hpp>
00055 #include <fei_Vector_Impl.hpp>
00056 #include <fei_Matrix_Impl.hpp>
00057 
00058 #ifdef HAVE_FEI_AZTECOO
00059 #include <fei_Aztec_LinSysCore.hpp>
00060 #endif
00061 #include <fei_Factory_Trilinos.hpp>
00062 
00063 #ifdef HAVE_FEI_FETI
00064 #include <FETI_DP_FiniteElementData.h>
00065 #endif
00066 
00067 #undef fei_file
00068 #define fei_file "test_Matrix.cpp"
00069 #include <fei_ErrMacros.hpp>
00070 
00071 int test_matrix_unit1()
00072 {
00073   std::vector<double> data1D;
00074   std::vector<const double*>data2D;
00075 
00076   int numRows = 2;
00077   int numCols = 3;
00078   double** values = new double*[numRows];
00079   int i, j;
00080   for(i=0; i<numRows; ++i) {
00081     values[i] = new double[numCols];
00082     for(j=0; j<numCols; ++j) {
00083       values[i][j] = i*1.0;
00084     }
00085   }
00086 
00087   fei::Matrix_core::copyTransposeToWorkArrays(numRows, numCols,
00088              values, data1D, data2D);
00089 
00090   for(i=0; i<numRows; ++i) {
00091     for(j=0; j<numCols; ++j) {
00092       if (std::abs(values[i][j] - (data2D[j])[i]) > 1.e-49) {
00093   ERReturn(-1);
00094       }
00095     }
00096     delete [] values[i];
00097   }
00098 
00099   delete [] values;
00100 
00101   return(0);
00102 }
00103 
00104 void test_Matrix_unit2(MPI_Comm comm, int numProcs, int localProc)
00105 {
00106   if (numProcs > 1) {
00107     return;
00108   }
00109 
00110   FEI_COUT << "testing fei::Matrix_Impl...";
00111 
00112   fei::SharedPtr<fei::VectorSpace> rowspace(new fei::VectorSpace(comm));
00113   fei::SharedPtr<fei::VectorSpace> colspace;
00114 
00115   int rowfield = 0, rowfieldsize = 1;
00116   int idType = 0;
00117   rowspace->defineFields(1, &rowfield, &rowfieldsize);
00118   rowspace->defineIDTypes(1, &idType);
00119 
00120   fei::SharedPtr<fei::MatrixGraph> mgraph(new fei::MatrixGraph_Impl2(rowspace, colspace));
00121 
00122   int patternID1 = mgraph->definePattern(2, idType, rowfield);
00123 
00124   fei::Pattern* rowpattern = mgraph->getPattern(patternID1);
00125 
00126   mgraph->initConnectivityBlock(0, 1, patternID1);
00127 
00128   std::vector<int> ids(2);
00129   ids[0] = 0; ids[1] = 1;
00130 
00131   int err = mgraph->initConnectivity(0, 0, &ids[0]);
00132   if (err) {
00133     FEI_OSTRINGSTREAM osstr;
00134     osstr << "test_Matrix_unit2, initConnectivity returned err="<<err;
00135     throw std::runtime_error(osstr.str());
00136   }
00137 
00138   err = mgraph->initComplete();
00139   if (err) {
00140     FEI_OSTRINGSTREAM osstr;
00141     osstr << "test_Matrix_unit2, initComplete returned err="<<err;
00142     throw std::runtime_error(osstr.str());
00143   }
00144 
00145   bool factory_created = false;
00146   fei::SharedPtr<fei::Factory> factory;
00147   try {
00148     factory = fei::create_fei_Factory(comm, "Trilinos");
00149     factory_created = true;
00150   }
00151   catch(...) {}
00152 
00153   if (!factory_created) {
00154     FEI_COUT << "failed to create Trilinos factory."<<FEI_ENDL;
00155     return;
00156   }
00157 
00158   fei::SharedPtr<fei::Matrix> feimat = factory->createMatrix(mgraph);
00159 
00160   int numrowindices = rowpattern->getNumIndices();
00161 
00162   std::vector<double> coefs(numrowindices*numrowindices, 1.0);
00163   std::vector<double*> coefs_2D(numrowindices);
00164   for(int i=0; i<numrowindices; ++i) {
00165     coefs_2D[i] = &(coefs[i*numrowindices]);
00166   }
00167 
00168   err = feimat->sumIn(0, 0, &coefs_2D[0]);
00169   if (err) {
00170     FEI_OSTRINGSTREAM osstr;
00171     osstr << "test_Matrix_unit2, feimat->sumIn returned err="<<err;
00172     throw std::runtime_error(osstr.str());
00173   }
00174 
00175   err = feimat->globalAssemble();
00176   if (err) {
00177     FEI_OSTRINGSTREAM osstr;
00178     osstr << "test_Matrix_unit2, feimat->globalAssemble returned err="<<err;
00179     throw std::runtime_error(osstr.str());
00180   }
00181 
00182   err = feimat->writeToFile("feimat2.mtx", false);
00183   if (err) {
00184     FEI_OSTRINGSTREAM osstr;
00185     osstr << "test_Matrix_unit2, feimat->writeToFile returned err="<<err;
00186     throw std::runtime_error(osstr.str());
00187   }
00188 
00189   fei::FillableMat feimat_ss;
00190   err = fei_test_utils::copy_feiMatrix_to_FillableMat(*feimat, feimat_ss);
00191   if (err) {
00192     FEI_OSTRINGSTREAM osstr;
00193     osstr << "test_Matrix_unit2, copy_feiMatrix_to_FillableMat returned err="<<err;
00194     throw std::runtime_error(osstr.str());
00195   }
00196 
00197   fei_test_utils::writeMatrix("feimat_ss2.mtx", feimat_ss);
00198 
00199   FEI_COUT << "ok"<<FEI_ENDL;
00200 }
00201 
00202 void test_Matrix_unit4(MPI_Comm comm, int numProcs, int localProc)
00203 {
00204   if (numProcs > 1) {
00205     return;
00206   }
00207 
00208   FEI_COUT << "testing fei::Matrix_Impl with FEI_BLOCK_DIAGONAL_ROW...";
00209 
00210   fei::SharedPtr<fei::VectorSpace> rowspace(new fei::VectorSpace(comm));
00211   fei::SharedPtr<fei::VectorSpace> colspace;
00212 
00213   int rowfield = 0, rowfieldsize = 2;
00214   int idType = 0;
00215   rowspace->defineFields(1, &rowfield, &rowfieldsize);
00216   rowspace->defineIDTypes(1, &idType);
00217 
00218   fei::SharedPtr<fei::MatrixGraph> mgraph(new fei::MatrixGraph_Impl2(rowspace, colspace));
00219 
00220   int patternID1 = mgraph->definePattern(2, idType, rowfield);
00221 
00222   fei::Pattern* rowpattern = mgraph->getPattern(patternID1);
00223 
00224   bool diagonal = true;
00225   mgraph->initConnectivityBlock(0, 1, patternID1, diagonal);
00226 
00227   std::vector<int> ids(2);
00228   ids[0] = 0; ids[1] = 1;
00229 
00230   int err = mgraph->initConnectivity(0, 0, &ids[0]);
00231   if (err) {
00232     FEI_OSTRINGSTREAM osstr;
00233     osstr << "test_Matrix_unit4, initConnectivity returned err="<<err;
00234     throw std::runtime_error(osstr.str());
00235   }
00236 
00237   err = mgraph->initComplete();
00238   if (err) {
00239     FEI_OSTRINGSTREAM osstr;
00240     osstr << "test_Matrix_unit4, initComplete returned err="<<err;
00241     throw std::runtime_error(osstr.str());
00242   }
00243 
00244   fei::SharedPtr<fei::Factory> factory;
00245   try {
00246     factory = fei::create_fei_Factory(comm, "Trilinos");
00247   }
00248   catch(...) {
00249     FEI_COUT << "Trilinos not available."<<FEI_ENDL;
00250     return;
00251   }
00252 
00253   fei::Param blktrue("BLOCK_MATRIX", true);
00254   fei::Param blkfalse("BLOCK_MATRIX", false);
00255 
00256   fei::ParameterSet paramset;
00257   paramset.add(blktrue);
00258   factory->parameters(paramset);
00259 
00260   fei::SharedPtr<fei::Matrix> feiblkmat = factory->createMatrix(mgraph);
00261 
00262   paramset.add(blkfalse);
00263   factory->parameters(paramset);
00264 
00265   fei::SharedPtr<fei::Matrix> feimat = factory->createMatrix(mgraph);
00266 
00267   int numrowindices = rowpattern->getNumIndices();
00268 
00269   std::vector<double> coefs(numrowindices*rowfieldsize*rowfieldsize, 1.0);
00270   std::vector<double*> coefs_2D(numrowindices*rowfieldsize);
00271   int offset = 0;
00272   for(int i=0; i<numrowindices*rowfieldsize; ++i) {
00273     coefs_2D[i] = &(coefs[offset]);
00274     offset += rowfieldsize;
00275   }
00276 
00277   err = feimat->sumIn(0, 0, &coefs_2D[0], FEI_BLOCK_DIAGONAL_ROW);
00278   if (err) {
00279     FEI_OSTRINGSTREAM osstr;
00280     osstr << "test_Matrix_unit4, feimat->sumIn returned err="<<err;
00281     throw std::runtime_error(osstr.str());
00282   }
00283 
00284   err = feiblkmat->sumIn(0, 0, &coefs_2D[0], FEI_BLOCK_DIAGONAL_ROW);
00285   if (err) {
00286     FEI_OSTRINGSTREAM osstr;
00287     osstr << "test_Matrix_unit4, feiblkmat->sumIn returned err="<<err;
00288     throw std::runtime_error(osstr.str());
00289   }
00290 
00291   err = feimat->globalAssemble();
00292   if (err) {
00293     FEI_OSTRINGSTREAM osstr;
00294     osstr << "test_Matrix_unit4, feimat->globalAssemble returned err="<<err;
00295     throw std::runtime_error(osstr.str());
00296   }
00297 
00298   err = feiblkmat->globalAssemble();
00299   if (err) {
00300     FEI_OSTRINGSTREAM osstr;
00301     osstr << "test_Matrix_unit4, feimat->globalAssemble returned err="<<err;
00302     throw std::runtime_error(osstr.str());
00303   }
00304 
00305   feimat->writeToFile("feimat_blkdiag.mtx");
00306   feiblkmat->writeToFile("feiblkmat_blkdiag.mtx");
00307 
00308   FEI_COUT << "ok"<<FEI_ENDL;
00309 }
00310 
00311 test_Matrix::test_Matrix(MPI_Comm comm)
00312   : tester(comm)
00313 {
00314 }
00315 
00316 test_Matrix::~test_Matrix()
00317 {
00318 }
00319 
00320 int test_Matrix::runtests()
00321 {
00322 
00323   //-------------------------------
00324   // Test a Factory_Trilinos matrix.
00325   if (localProc_==0) FEI_COUT << "Testing Factory_Trilinos fei::Matrix..." << FEI_ENDL;
00326   fei::SharedPtr<fei::Factory> factory_trilinos(new Factory_Trilinos(comm_));
00327 
00328   fei::SharedPtr<fei::Matrix> mat = create_matrix(factory_trilinos);
00329 
00330   matrix_test1(mat);
00331 
00332   if (localProc_==0) FEI_COUT << FEI_ENDL;
00333 
00334 #ifdef HAVE_FEI_AZTECOO
00335   //-------------------------------
00336   // Test a matrix from a "snl_fei::Factory", which was created by the
00337   // test-util function create_fei_Factory for library "Aztec", which causes
00338   // an "old" LinearSystemCore object to be used as the underlying assembly layer.
00339   // (It still ends up using Aztec if you create a solver and solve the linear-
00340   // system, but it assembles data into non-Trilinos matrix/vector data structures.)
00341   if (localProc_==0) FEI_COUT << "Testing 'old' factory_aztec fei::Matrix..." << FEI_ENDL;
00342   fei::SharedPtr<fei::Factory> factory_aztec =
00343     fei::create_fei_Factory(comm_, "Aztec");
00344 
00345   fei::SharedPtr<fei::Matrix> mat1 = create_matrix(factory_aztec);
00346 
00347   matrix_test1(mat1);
00348 #endif
00349 
00350   return(0);
00351 }
00352 
00353 fei::SharedPtr<fei::Matrix>
00354 test_Matrix::create_matrix(fei::SharedPtr<fei::Factory> factory)
00355 {
00356   testData test_data(localProc_, numProcs_);
00357 
00358   fei::SharedPtr<fei::VectorSpace> vspace =
00359     test_VectorSpace::create_VectorSpace(comm_, &test_data, localProc_, numProcs_,
00360            false, false, (const char*)0, factory);
00361   int err = vspace->initComplete();
00362   if (err != 0) {
00363     FEI_COUT << "ERROR, failed to create valid fei::VectorSpace." << FEI_ENDL;
00364     throw std::runtime_error("test_Vector::vector_test1: ERROR, failed to create valid fei::VectorSpace.");
00365   }
00366 
00367   fei::SharedPtr<fei::MatrixGraph> mgraph =
00368     factory->createMatrixGraph(vspace, vspace, NULL);
00369 
00370   std::vector<int>& fieldIDs = test_data.fieldIDs;
00371   std::vector<int>& idTypes = test_data.idTypes;
00372   std::vector<int>& ids = test_data.ids;
00373 
00374   int numIDs = ids.size();
00375   int fieldID = fieldIDs[0];
00376   int idType = idTypes[0];
00377 
00378   int patternID = mgraph->definePattern(numIDs, idType, fieldID);
00379 
00380   mgraph->initConnectivityBlock(0, 1, patternID);
00381 
00382   mgraph->initConnectivity(0, 0, &ids[0]);
00383 
00384   mgraph->initComplete();
00385 
00386   fei::SharedPtr<fei::Matrix> matrix = factory->createMatrix(mgraph);
00387   return(matrix);
00388 }
00389 
00390 void test_Matrix::matrix_test1(fei::SharedPtr<fei::Matrix> mat)
00391 {
00392   if (localProc_==0)
00393     FEI_COUT << "  matrix_test1: testing fei::Matrix with type '"
00394      << mat->typeName() << "':"<<FEI_ENDL;
00395 
00396   fei::SharedPtr<fei::MatrixGraph> mgraph = mat->getMatrixGraph();
00397 
00398   fei::SharedPtr<fei::VectorSpace> rspace = mgraph->getRowSpace();
00399 
00400   if (localProc_==0)
00401     FEI_COUT << "   testing get{Global/Local}NumRows,getRowLength...";
00402 
00403   int mglobalrows = mat->getGlobalNumRows();
00404   int vglobaleqns = rspace->getGlobalNumIndices();
00405 
00406   if (mglobalrows != vglobaleqns) {
00407     throw std::runtime_error("mat reports different num rows than vector-space eqns");
00408   }
00409 
00410   std::vector<int> global_offsets;
00411   rspace->getGlobalIndexOffsets(global_offsets);
00412 
00413   int my_num_rows = mat->getLocalNumRows();
00414   if (my_num_rows != global_offsets[localProc_+1]-global_offsets[localProc_]) {
00415     throw std::runtime_error("num-local-rows mis-match between mat and vector-space");
00416   }
00417 
00418   int i, my_first_row = global_offsets[localProc_];
00419   std::vector<int> row_lengths(my_num_rows);
00420 
00421   for(i=0; i<my_num_rows; ++i) {
00422     int errcode = mat->getRowLength(i+my_first_row, row_lengths[i]);
00423     if (errcode != 0) {
00424       throw std::runtime_error("nonzero errcode from mat->getRowLength");
00425     }
00426   }
00427 
00428   if (localProc_==0) FEI_COUT << "ok" << FEI_ENDL;
00429 }
00430 
00431 int test_Matrix::serialtest1()
00432 {
00433   fei::SharedPtr<testData> testdata(new testData(localProc_, numProcs_));
00434   std::vector<int>& fieldIDs = testdata->fieldIDs;
00435   std::vector<int>& fieldSizes = testdata->fieldSizes;
00436   std::vector<int>& idTypes = testdata->idTypes;
00437   std::vector<int>& ids = testdata->ids;
00438 
00439   fei::SharedPtr<fei::VectorSpace> vspc(new fei::VectorSpace(comm_, "sU_Mat"));
00440 
00441   vspc->defineFields(fieldIDs.size(),
00442          &fieldIDs[0],
00443          &fieldSizes[0]);
00444 
00445   vspc->defineIDTypes(idTypes.size(), &idTypes[0]);
00446 
00447   fei::SharedPtr<fei::MatrixGraph>
00448     matgraph(new fei::MatrixGraph_Impl2(vspc, vspc, "sU_Mat"));
00449 
00450   int numIDs = ids.size();
00451   int fieldID = fieldIDs[0];
00452   int idType = idTypes[0];
00453 
00454   int patternID = matgraph->definePattern(numIDs, idType, fieldID);
00455 
00456   CHK_ERR( matgraph->initConnectivityBlock(0, 1, patternID) );
00457 
00458   CHK_ERR( matgraph->initConnectivity(0, 0, &ids[0]) );
00459 
00460   CHK_ERR( matgraph->initComplete() );
00461 
00462   fei::SharedPtr<fei::FillableMat> ssmat(new fei::FillableMat), ssmatT(new fei::FillableMat);
00463   int localsize = matgraph->getRowSpace()->getNumIndices_Owned();
00464   fei::SharedPtr<fei::Matrix> matrix(new fei::Matrix_Impl<fei::FillableMat>(ssmat, matgraph, localsize));
00465 
00466   fei::SharedPtr<fei::Matrix> matrixT(new fei::Matrix_Impl<fei::FillableMat>(ssmatT, matgraph, localsize));
00467 
00468   std::vector<int> indices(numIDs);
00469   CHK_ERR( matgraph->getConnectivityIndices(0, 0, numIDs, &indices[0], numIDs) );
00470 
00471   std::vector<double> data1(numIDs*numIDs);
00472   std::vector<double*> data2d(numIDs);
00473 
00474   int i;
00475   for(i=0; i<numIDs; ++i) {
00476     data2d[i] = &(data1[i*numIDs]);
00477   }
00478 
00479   for(i=0; i<numIDs*numIDs; ++i) {
00480     data1[i] = 1.0*i;
00481   }
00482 
00483   CHK_ERR( matrix->sumIn(numIDs, &indices[0], numIDs, &indices[0],
00484        &data2d[0], 0) );
00485 
00486   CHK_ERR( matrix->sumIn(0, 0, &data2d[0], 0) );
00487 
00488   CHK_ERR( matrixT->sumIn(numIDs, &indices[0],
00489        numIDs, &indices[0], &data2d[0], 3) );
00490 
00491   CHK_ERR( matrixT->sumIn(0, 0, &data2d[0], 3) );
00492 
00493   if (*ssmat != *ssmatT) {
00494     ERReturn(-1);
00495   }
00496 
00497   return(0);
00498 }
00499 
00500 int test_Matrix::serialtest2()
00501 {
00502   testData* testdata = new testData(localProc_, numProcs_);
00503   std::vector<int>& idTypes = testdata->idTypes;
00504   std::vector<int>& ids = testdata->ids;
00505 
00506   fei::SharedPtr<fei::VectorSpace> vspc(new fei::VectorSpace(comm_, "sU_Mat"));
00507 
00508   vspc->defineIDTypes(idTypes.size(), &idTypes[0]);
00509 
00510   fei::SharedPtr<fei::MatrixGraph> matgraph(new fei::MatrixGraph_Impl2(vspc, vspc, "sU_Mat"));
00511 
00512   int numIDs = ids.size();
00513   int idType = idTypes[0];
00514 
00515   int patternID = matgraph->definePattern(numIDs, idType);
00516 
00517   CHK_ERR( matgraph->initConnectivityBlock(0, 1, patternID) );
00518 
00519   CHK_ERR( matgraph->initConnectivity(0, 0, &ids[0]) );
00520 
00521   CHK_ERR( matgraph->initComplete() );
00522 
00523   fei::SharedPtr<fei::FillableMat> ssmat(new fei::FillableMat), ssmatT(new fei::FillableMat);
00524   int localsize = matgraph->getRowSpace()->getNumIndices_Owned();
00525   fei::Matrix* matrix = new fei::Matrix_Impl<fei::FillableMat>(ssmat, matgraph, localsize);
00526 
00527   fei::Matrix* matrixT = new fei::Matrix_Impl<fei::FillableMat>(ssmatT, matgraph, localsize);
00528 
00529   std::vector<int> indices(numIDs);
00530   CHK_ERR( matgraph->getConnectivityIndices(0, 0, numIDs, &indices[0], numIDs) );
00531 
00532   std::vector<double> data1(numIDs*numIDs);
00533   std::vector<double*> data2d(numIDs);
00534 
00535   int i;
00536   for(i=0; i<numIDs; ++i) {
00537     data2d[i] = &(data1[i*numIDs]);
00538   }
00539 
00540   for(i=0; i<numIDs*numIDs; ++i) {
00541     data1[i] = 1.0*i;
00542   }
00543 
00544   CHK_ERR( matrix->sumIn(numIDs, &indices[0],
00545        numIDs, &indices[0], &data2d[0], 0) );
00546 
00547   CHK_ERR( matrixT->sumIn(numIDs, &indices[0],
00548        numIDs, &indices[0], &data2d[0], 3) );
00549 
00550   if (*ssmat != *ssmatT) {
00551     ERReturn(-1);
00552   }
00553 
00554   delete matrix;
00555   delete matrixT;
00556   delete testdata;
00557 
00558   return(0);
00559 }
00560 
00561 int test_Matrix::serialtest3()
00562 {
00563   testData* testdata = new testData(localProc_, numProcs_);
00564   std::vector<int>& fieldIDs = testdata->fieldIDs;
00565   std::vector<int>& fieldSizes = testdata->fieldSizes;
00566   std::vector<int>& idTypes = testdata->idTypes;
00567   std::vector<int>& ids = testdata->ids;
00568 
00569   fei::SharedPtr<fei::VectorSpace> vspc(new fei::VectorSpace(comm_, "sU_Mat3"));
00570 
00571   vspc->defineFields(fieldIDs.size(), &fieldIDs[0], &fieldSizes[0]);
00572 
00573   vspc->defineIDTypes(idTypes.size(), &idTypes[0]);
00574 
00575   fei::SharedPtr<fei::MatrixGraph>
00576     matgraph(new fei::MatrixGraph_Impl2(vspc, vspc, "sU_Mat3"));
00577 
00578   int numIDs = ids.size();
00579   int fieldID = fieldIDs[0];
00580   int idType = idTypes[0];
00581 
00582   int patternID = matgraph->definePattern(numIDs, idType, fieldID);
00583 
00584   CHK_ERR( matgraph->initConnectivityBlock(0, 1, patternID) );
00585 
00586   CHK_ERR( matgraph->initConnectivity(0, 0, &ids[0]) );
00587 
00588   //set up a slave constraint that defines id 2, field 0 to be equal to
00589   //id 1, field 0.
00590   int offsetOfSlave = 1;
00591   int offsetIntoSlaveField = 0;
00592   std::vector<double> weights(2);
00593   weights[0] = 1.0;
00594   weights[1] = -1.0;
00595   double rhsValue = 0.0;
00596   std::vector<int> cr_idtypes(2, idTypes[0]);
00597   std::vector<int> cr_fieldIDs(2, fieldIDs[0]);
00598 
00599   CHK_ERR( matgraph->initSlaveConstraint(2, //numIDs
00600           &cr_idtypes[0],
00601           &ids[1],
00602           &cr_fieldIDs[0],
00603           offsetOfSlave,
00604           offsetIntoSlaveField,
00605           &weights[0],
00606           rhsValue) );
00607 
00608   CHK_ERR( matgraph->initComplete() );
00609 
00610   fei::SharedPtr<fei::FillableMat> ssmat(new fei::FillableMat);
00611   int localsize = matgraph->getRowSpace()->getNumIndices_Owned();
00612   localsize -= 1;//subtract the slave
00613   fei::Matrix* matrix = new fei::Matrix_Impl<fei::FillableMat>(ssmat, matgraph, localsize);
00614 
00615   if (matrix == NULL) {
00616     ERReturn(-1);
00617   }
00618 
00619   std::vector<int> indices(numIDs);
00620   CHK_ERR( matgraph->getConnectivityIndices(0, 0, numIDs,
00621              &indices[0], numIDs) );
00622 
00623   std::vector<double> data1(numIDs*numIDs);
00624   std::vector<double*> data2d(numIDs);
00625 
00626   int i;
00627   for(i=0; i<numIDs; ++i) {
00628     data2d[i] = &(data1[i*numIDs]);
00629   }
00630 
00631   for(i=0; i<numIDs*numIDs; ++i) {
00632     data1[i] = 1.0*i;
00633   }
00634 
00635   CHK_ERR( matrix->sumIn(numIDs, &indices[0],
00636        numIDs, &indices[0], &data2d[0], 0) );
00637 
00638   CHK_ERR( matrix->sumIn(0, 0, &data2d[0], 0) );
00639 
00640   delete matrix;
00641   delete testdata;
00642 
00643   return(0);
00644 }
00645 
00646 int test_Matrix::test1()
00647 {
00648   return(0);
00649 }
00650 
00651 int test_Matrix::test2()
00652 {
00653  return(0);
00654 }
00655 
00656 int test_Matrix::test3()
00657 {
00658 #ifdef HAVE_FEI_FETI
00659   testData* testdata = new testData(localProc_, numProcs_);
00660   std::vector<int>& idTypes = testdata->idTypes;
00661   std::vector<int>& ids = testdata->ids;
00662 
00663   fei::SharedPtr<FiniteElementData> fedata(new FETI_DP_FiniteElementData(comm_));
00664 
00665   std::string paramstr("debugOutput .");
00666   char* param = const_cast<char*>(paramstr.c_str());
00667 
00668   CHK_ERR( fedata->parameters(1, &param) );
00669 
00670   fei::SharedPtr<fei::Factory> factory(new snl_fei::Factory(fedata, idTypes[0]));
00671 
00672   fei::SharedPtr<fei::VectorSpace> vectorSpacePtr =
00673     test_VectorSpace::create_VectorSpace(comm_,
00674            testdata, localProc_, numProcs_,
00675            false, false, "U_FEMat", factory);
00676 
00677   fei::SharedPtr<fei::MatrixGraph> matrixGraphPtr =
00678     test_MatrixGraph::create_MatrixGraph(testdata, localProc_, numProcs_,
00679            false, false, "U_FEMat", vectorSpacePtr,
00680            factory);
00681 
00682   CHK_ERR( matrixGraphPtr->initComplete() );
00683 
00684   fei::SharedPtr<fei::Vector> vec_fed = factory->createVector(vectorSpacePtr);
00685 
00686   fei::SharedPtr<fei::Matrix> mat_fed = factory->createMatrix(matrixGraphPtr);
00687 
00688   fei::Matrix_Impl<FiniteElementData>* smat2 = 
00689     dynamic_cast<fei::Matrix_Impl<FiniteElementData>*>(mat_fed.get());
00690   if (smat2 == NULL) {
00691     ERReturn(-1);
00692   }
00693 
00694   int blockID=0;
00695   int numIndices = matrixGraphPtr->getConnectivityNumIndices(blockID);
00696 
00697   std::vector<int> indicesArray(numIndices);
00698   int* indicesPtr = &indicesArray[0];
00699 
00700   int checkNumIndices = 0;
00701   CHK_ERR( matrixGraphPtr->getConnectivityIndices(blockID, 0,
00702                numIndices, indicesPtr,
00703                checkNumIndices) );
00704 
00705   std::vector<double> data(ids.size(), 1.0);
00706   double* dptr = &data[0];
00707   std::vector<double*> coefPtrs(ids.size(), dptr);
00708 
00709   CHK_ERR( mat_fed->sumIn(blockID, 0, &coefPtrs[0]) );
00710 
00711   CHK_ERR( vec_fed->sumIn(blockID, 0, &data[0]) );
00712 
00713   CHK_ERR( mat_fed->gatherFromOverlap() );
00714 
00715   CHK_ERR( fedata->loadComplete() );
00716 
00717 
00718   delete testdata;
00719 
00720   MPI_Barrier(comm_);
00721 
00722 #endif  //HAVE_FEI_FETI
00723 
00724   return(0);
00725 }
00726 
00727 int test_Matrix::test4()
00728 {
00729   return(0);
00730 }
00731 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends