Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Kokkos_Mkl_MatrixDescriptor.cpp
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //          Kokkos: Node API and Parallel Node Kernels
00005 //              Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the 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 Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #include "Kokkos_Mkl_MatrixDescriptor.hpp"
00043 
00044 
00045 namespace Kokkos {
00046   namespace Mkl {
00047     void MatrixDescriptor::print (Teuchos::FancyOStream& out) const
00048     {
00049       using std::endl;
00050       TEUCHOS_TEST_FOR_EXCEPTION(descr_.size() < 6, std::invalid_argument,
00051         "MKL documentation requires that the sparse matrix descriptor array "
00052         "have at least 6 elements.  MKL only uses the first 4 currently.");
00053 
00054       const char s0 = descr_[0];
00055       std::string structure;
00056       if (s0 == 'G') {
00057         structure = "General";
00058       } else if (s0 == 'S') {
00059         structure = "Symmetric";
00060       } else if (s0 == 'H') {
00061         structure = "Hermitian";
00062       } else if (s0 == 'T') {
00063         structure = "Triangular";
00064       } else if (s0 == 'A') {
00065         structure = "Antisymmetric";
00066       } else if (s0 == 'D') {
00067         structure = "Diagonal";
00068       } else {
00069         TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00070           "First element of the matrix descriptor array is invalid.  "
00071           "Valid values are 'G', 'S', 'H', 'T', 'A', and 'D', but you "
00072           "provided '" << s0 << "'.");
00073       }
00074 
00075       // Symmetric, Hermitian, and Triangular matrices all use
00076       // triangular storage, in the sense that MKL ignores all entries
00077       // in the lower resp. upper triangle of the matrix in these
00078       // cases.
00079       const char s1 = descr_[1];
00080       std::string upperLower;
00081       if (s1 == 'L') {
00082         upperLower = "Lower";
00083       } else if (s1 == 'U') {
00084         upperLower = "Upper";
00085       } else {
00086         upperLower = "Neither";
00087       }
00088 
00089       const char s2 = descr_[2];
00090       const bool unitDiag = (s2 == 'U');
00091 
00092       // 'C' means zero-based indexing (as in C); 'F' means one-based
00093       // indexing (as in Fortran).
00094       const char s3 = descr_[3];
00095       int indexBase = 0;
00096       if (s3 == 'C') {
00097         indexBase = 0;
00098       } else if (s3 == 'F') {
00099         indexBase = 1;
00100       } else {
00101         TEUCHOS_TEST_FOR_EXCEPTION(
00102           s3 != 'C' && s3 != 'F',
00103           std::invalid_argument,
00104           "Fourth element of matrix descriptor has invalid index base "
00105           "specification '" << s3 << "'.  Valid values are 'C' (for zero-based, "
00106           "C-style indices) and 'F' (for one-based, Fortran-style indices).");
00107       }
00108 
00109       // YAML (Yet Another Markup Language) output.
00110       // We defer all output until end of routine, for exception safety.
00111       out << "MKL sparse matrix descriptor:" << endl;
00112       Teuchos::OSTab tab (Teuchos::rcpFromRef (out));
00113       out << "Structure: \"" << structure << "\"" << endl
00114           << "UpperLower: \"" << upperLower << "\"" << endl
00115           << "UnitDiagonal: " << (unitDiag ? "true" : "false") << endl
00116           << "IndexBase: " << indexBase << endl;
00117     }
00118 
00119     MatrixDescriptor::
00120     MatrixDescriptor (const bool symmetric,
00121                       const bool hermitian,
00122                       const bool triangular,
00123                       const bool skewSymmetric,
00124                       const bool diagonal,
00125                       const Teuchos::EUplo uplo,
00126                       const Teuchos::EDiag diag,
00127                       const int indexBase) :
00128       descr_ ("XXXXXX") // initialize with correct length but invalid values
00129     {
00130       fill (descr_, symmetric, hermitian, triangular, skewSymmetric,
00131             diagonal, uplo, diag, indexBase);
00132     }
00133 
00134     MatrixDescriptor::MatrixDescriptor () :
00135       descr_ ("XXXXXX") // initialize with correct length but invalid values
00136     {
00137       const bool symmetric = false;
00138       const bool hermitian = false;
00139       const bool triangular = false;
00140       const bool skewSymmetric = false;
00141       const bool diagonal = false;
00142       const Teuchos::EUplo uplo = Teuchos::UNDEF_TRI;
00143       const Teuchos::EDiag diag = Teuchos::NON_UNIT_DIAG;
00144       // We use one-based indices by default, so that MKL's mat-vec
00145       // routines that take multiple vectors accept the vectors in
00146       // column-major order.  With zero-based indices, those routines
00147       // require the vectors in row-major order.
00148       const int indexBase = 1;
00149       fill (descr_, symmetric, hermitian, triangular, skewSymmetric,
00150             diagonal, uplo, diag, indexBase);
00151     }
00152 
00153     void
00154     MatrixDescriptor::fill (std::string& descr,
00155                             const bool symmetric,
00156                             const bool hermitian,
00157                             const bool triangular,
00158                             const bool skewSymmetric,
00159                             const bool diagonal,
00160                             const Teuchos::EUplo uplo,
00161                             const Teuchos::EDiag diag,
00162                             const int indexBase)
00163     {
00164       // MKL docs require 6 elements, of which only the first 4 are
00165       // used currently.
00166       if (descr.size() < 6) {
00167         descr.resize (6);
00168         // Just put something here that's different, as a marker.
00169         // MKL doesn't use it, so we can put whatever we want.
00170         descr[4] = 'X';
00171         descr[5] = 'X';
00172       }
00173 
00174       char s0 = 'G';
00175       // Diagonal test goes first, since a diagonal matrix may also be
00176       // symmetric, Hermitian, or skew-symmetric, depending on the
00177       // diagonal's values.  The user is responsible for verifying
00178       // symmetry, Hermitian-ness, or skew-symmetry of a diagonal
00179       // matrix.
00180       if (diagonal) {
00181         s0 = 'D';
00182       } else if (symmetric) {
00183         s0 = 'S';
00184       } else if (hermitian) {
00185         s0 = 'H';
00186       } else if (uplo != Teuchos::UNDEF_TRI) {
00187         s0 = 'T';
00188       } else if (skewSymmetric) {
00189         s0 = 'A'; // A for antisymmetric
00190       } else {
00191         s0 = 'G'; // G for general
00192       }
00193 
00194       // If the matrix is neither upper nor lower triangular, we put
00195       // 'N' (for "Neither").  Note that MKL's routines dealing with
00196       // triangular, (skew-)symmetric, or Hermitian matrices only read
00197       // the elements of the matrix in the lower resp. upper triangle,
00198       // even if storage includes both lower and upper triangles.
00199       char s1 = 'N';
00200       if (uplo == Teuchos::LOWER_TRI) {
00201         s1 = 'L';
00202       } else if (uplo == Teuchos::UPPER_TRI) {
00203         s1 = 'U';
00204       } else {
00205         s1 = 'N';
00206       }
00207 
00208       // Whether or not the matrix has unit diagonal.
00209       char s2 = 'N';
00210       if (diag == Teuchos::NON_UNIT_DIAG) {
00211         s2 = 'N';
00212       } else if (diag == Teuchos::UNIT_DIAG) {
00213         s2 = 'U';
00214       } else {
00215         TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00216           "Invalid Teuchos::EDiag enum value " << diag << ".  Valid values are "
00217           "Teuchos::NON_UNIT_DIAG and Teuchos::UNIT_DIAG.");
00218       }
00219 
00220       // 'C' means zero-based indexing (as in C); 'F' means one-based
00221       // indexing (as in Fortran).
00222       char s3 = 'C';
00223       if (indexBase == 0) {
00224         s3 = 'C';
00225       } else if (indexBase == 1) {
00226         s3 = 'F';
00227       } else {
00228         TEUCHOS_TEST_FOR_EXCEPTION(indexBase != 0 && indexBase != 1,
00229           std::invalid_argument, "Invalid index base " << indexBase
00230           << ".  Valid values are 0 (for zero-based indexing) and 1 "
00231           "(for one-based indexing).");
00232       }
00233 
00234       // Now that we've validated the input, commit changes.
00235       descr[0] = s0;
00236       descr[1] = s1;
00237       descr[2] = s2;
00238       descr[3] = s3;
00239     }
00240 
00241     int
00242     MatrixDescriptor::getIndexBase () const
00243     {
00244       if (descr_[3] == 'C') {
00245         return 0;
00246       } else if (descr_[3] == 'F') {
00247         return 1;
00248       } else {
00249         TEUCHOS_TEST_FOR_EXCEPTION(
00250           true, std::invalid_argument, "Invalid MKL matrix descriptor: fourth "
00251           "entry of the character array must be either 'C' or 'F', but is '"
00252           << descr_[3] << "' instead.");
00253       }
00254     }
00255 
00256     Teuchos::EDiag
00257     MatrixDescriptor::getDiag () const
00258     {
00259       const char s2 = descr_[2];
00260       if (s2 == 'U') {
00261         return Teuchos::UNIT_DIAG;
00262       } else {
00263         return Teuchos::NON_UNIT_DIAG;
00264       }
00265     }
00266 
00267     Teuchos::EUplo
00268     MatrixDescriptor::getUplo () const
00269     {
00270       const char s1 = descr_[1];
00271       if (s1 == 'L') {
00272         return Teuchos::LOWER_TRI;
00273       } else if (s1 == 'U') {
00274         return Teuchos::UPPER_TRI;
00275       } else {
00276         return Teuchos::UNDEF_TRI;
00277       }
00278     }
00279 
00280   } // namespace Mkl
00281 } // namespace Kokkos
00282 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends