Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Kokkos_Raw_SparseTriangularSolve_def.hpp
Go to the documentation of this file.
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 #ifndef __Kokkos_Raw_SparseTriangularSolve_def_hpp
00043 #define __Kokkos_Raw_SparseTriangularSolve_def_hpp
00044 
00050 
00051 namespace Kokkos {
00052 namespace Raw {
00053 
00054 template<class Ordinal,
00055          class MatrixScalar,
00056          class DomainScalar,
00057          class RangeScalar>
00058 void
00059 lowerTriSolveCscColMajor (
00060   const Ordinal numRows,
00061   const Ordinal numCols,
00062   const Ordinal numVecs,
00063   RangeScalar* const X,
00064   const Ordinal colStrideX,
00065   const  size_t* const ptr,
00066   const Ordinal* const ind,
00067   const MatrixScalar* const val,
00068   const DomainScalar* const Y,
00069   const Ordinal colStrideY)
00070 {
00071   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00072 
00073   for (Ordinal j = 0; j < numVecs; ++j) {
00074     for (Ordinal i = 0; i < numRows; ++i) {
00075       X[i + j*colStrideX] = Y[i + j*colStrideY];
00076     }
00077   }
00078 
00079   for (Ordinal c = 0; c < numCols; ++c) {
00080     MatrixScalar A_cc = STS::zero ();
00081     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00082       const Ordinal r = ind[k];
00083       MatrixScalar A_rc = val[k];
00084       if (r == c) {
00085         A_cc = A_cc + A_rc;
00086       } else {
00087         for (Ordinal j = 0; j < numVecs; ++j) {
00088           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00089         }
00090       }
00091       for (Ordinal j = 0; j < numVecs; ++j) {
00092         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
00093       }
00094     }
00095   }
00096 }
00097 
00098 template<class Ordinal,
00099          class MatrixScalar,
00100          class DomainScalar,
00101          class RangeScalar>
00102 void
00103 lowerTriSolveCsrColMajor (
00104   const Ordinal numRows,
00105   const Ordinal numCols,
00106   const Ordinal numVecs,
00107   RangeScalar* const X,
00108   const Ordinal colStrideX,
00109   const  size_t* const ptr,
00110   const Ordinal* const ind,
00111   const MatrixScalar* const val,
00112   const DomainScalar* const Y,
00113   const Ordinal colStrideY)
00114 {
00115   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00116 
00117   for (Ordinal r = 0; r < numRows; ++r) {
00118     for (Ordinal j = 0; j < numVecs; ++j) {
00119       X[r + j*colStrideX] = Y[r + j*colStrideY];
00120     }
00121     // We assume the diagonal entry is first in the row.
00122     const MatrixScalar A_rr = val[ptr[r]];
00123     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
00124       const MatrixScalar A_rc = val[k];
00125       const Ordinal c = ind[k];
00126       for (Ordinal j = 0; j < numVecs; ++j) {
00127         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00128       }
00129     }
00130     for (Ordinal j = 0; j < numVecs; ++j) {
00131       X[r + j*colStrideX] = X[r + j*colStrideX] / A_rr;
00132     }
00133   }
00134 }
00135 
00136 template<class Ordinal,
00137          class MatrixScalar,
00138          class DomainScalar,
00139          class RangeScalar>
00140 void
00141 lowerTriSolveCscRowMajor (
00142   const Ordinal numRows,
00143   const Ordinal numCols,
00144   const Ordinal numVecs,
00145   RangeScalar* const X,
00146   const Ordinal rowStrideX,
00147   const  size_t* const ptr,
00148   const Ordinal* const ind,
00149   const MatrixScalar* const val,
00150   const DomainScalar* const Y,
00151   const Ordinal rowStrideY)
00152 {
00153   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00154 
00155   for (Ordinal i = 0; i < numRows; ++i) {
00156     for (Ordinal j = 0; j < numVecs; ++j) {
00157       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
00158     }
00159   }
00160 
00161   for (Ordinal c = 0; c < numCols; ++c) {
00162     MatrixScalar A_cc = STS::zero ();
00163     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00164       const Ordinal r = ind[k];
00165       MatrixScalar A_rc = val[k];
00166       if (r == c) {
00167         A_cc = A_cc + A_rc;
00168       } else {
00169         for (Ordinal j = 0; j < numVecs; ++j) {
00170           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00171         }
00172       }
00173       for (Ordinal j = 0; j < numVecs; ++j) {
00174         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
00175       }
00176     }
00177   }
00178 }
00179 
00180 template<class Ordinal,
00181          class MatrixScalar,
00182          class DomainScalar,
00183          class RangeScalar>
00184 void
00185 lowerTriSolveCsrRowMajor (
00186   const Ordinal numRows,
00187   const Ordinal numCols,
00188   const Ordinal numVecs,
00189   RangeScalar* const X,
00190   const Ordinal rowStrideX,
00191   const  size_t* const ptr,
00192   const Ordinal* const ind,
00193   const MatrixScalar* const val,
00194   const DomainScalar* const Y,
00195   const Ordinal rowStrideY)
00196 {
00197   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00198 
00199   for (Ordinal r = 0; r < numRows; ++r) {
00200     for (Ordinal j = 0; j < numVecs; ++j) {
00201       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
00202     }
00203     // We assume the diagonal entry is first in the row.
00204     const MatrixScalar A_rr = val[ptr[r]];
00205     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
00206       const MatrixScalar A_rc = val[k];
00207       const Ordinal c = ind[k];
00208       for (Ordinal j = 0; j < numVecs; ++j) {
00209         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00210       }
00211     }
00212     for (Ordinal j = 0; j < numVecs; ++j) {
00213       X[r*rowStrideX + j] = X[r*rowStrideX + j] / A_rr;
00214     }
00215   }
00216 }
00217 
00218 template<class Ordinal,
00219          class MatrixScalar,
00220          class DomainScalar,
00221          class RangeScalar>
00222 void
00223 lowerTriSolveCscColMajorUnitDiag (
00224   const Ordinal numRows,
00225   const Ordinal numCols,
00226   const Ordinal numVecs,
00227   RangeScalar* const X,
00228   const Ordinal colStrideX,
00229   const  size_t* const ptr,
00230   const Ordinal* const ind,
00231   const MatrixScalar* const val,
00232   const DomainScalar* const Y,
00233   const Ordinal colStrideY)
00234 {
00235   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00236 
00237   for (Ordinal j = 0; j < numVecs; ++j) {
00238     for (Ordinal i = 0; i < numRows; ++i) {
00239       X[i + j*colStrideX] = Y[i + j*colStrideY];
00240     }
00241   }
00242 
00243   for (Ordinal c = 0; c < numCols; ++c) {
00244     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00245       const Ordinal r = ind[k];
00246       MatrixScalar A_rc = val[k];
00247       for (Ordinal j = 0; j < numVecs; ++j) {
00248         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00249       }
00250     }
00251   }
00252 }
00253 
00254 template<class Ordinal,
00255          class MatrixScalar,
00256          class DomainScalar,
00257          class RangeScalar>
00258 void
00259 lowerTriSolveCsrColMajorUnitDiag (
00260   const Ordinal numRows,
00261   const Ordinal numCols,
00262   const Ordinal numVecs,
00263   RangeScalar* const X,
00264   const Ordinal colStrideX,
00265   const  size_t* const ptr,
00266   const Ordinal* const ind,
00267   const MatrixScalar* const val,
00268   const DomainScalar* const Y,
00269   const Ordinal colStrideY)
00270 {
00271   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00272 
00273   for (Ordinal r = 0; r < numRows; ++r) {
00274     for (Ordinal j = 0; j < numVecs; ++j) {
00275       X[r + j*colStrideX] = Y[r + j*colStrideY];
00276     }
00277     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
00278       const MatrixScalar A_rc = val[k];
00279       const Ordinal c = ind[k];
00280       for (Ordinal j = 0; j < numVecs; ++j) {
00281         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00282       }
00283     }
00284   }
00285 }
00286 
00287 template<class Ordinal,
00288          class MatrixScalar,
00289          class DomainScalar,
00290          class RangeScalar>
00291 void
00292 lowerTriSolveCscRowMajorUnitDiag (
00293   const Ordinal numRows,
00294   const Ordinal numCols,
00295   const Ordinal numVecs,
00296   RangeScalar* const X,
00297   const Ordinal rowStrideX,
00298   const  size_t* const ptr,
00299   const Ordinal* const ind,
00300   const MatrixScalar* const val,
00301   const DomainScalar* const Y,
00302   const Ordinal rowStrideY)
00303 {
00304   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00305 
00306   for (Ordinal i = 0; i < numRows; ++i) {
00307     for (Ordinal j = 0; j < numVecs; ++j) {
00308       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
00309     }
00310   }
00311 
00312   for (Ordinal c = 0; c < numCols; ++c) {
00313     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00314       const Ordinal r = ind[k];
00315       MatrixScalar A_rc = val[k];
00316       for (Ordinal j = 0; j < numVecs; ++j) {
00317         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00318       }
00319     }
00320   }
00321 }
00322 
00323 template<class Ordinal,
00324          class MatrixScalar,
00325          class DomainScalar,
00326          class RangeScalar>
00327 void
00328 lowerTriSolveCsrRowMajorUnitDiag (
00329   const Ordinal numRows,
00330   const Ordinal numCols,
00331   const Ordinal numVecs,
00332   RangeScalar* const X,
00333   const Ordinal rowStrideX,
00334   const  size_t* const ptr,
00335   const Ordinal* const ind,
00336   const MatrixScalar* const val,
00337   const DomainScalar* const Y,
00338   const Ordinal rowStrideY)
00339 {
00340   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00341 
00342   for (Ordinal r = 0; r < numRows; ++r) {
00343     for (Ordinal j = 0; j < numVecs; ++j) {
00344       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
00345     }
00346     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
00347       const MatrixScalar A_rc = val[k];
00348       const Ordinal c = ind[k];
00349       for (Ordinal j = 0; j < numVecs; ++j) {
00350         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00351       }
00352     }
00353   }
00354 }
00355 
00356 template<class Ordinal,
00357          class MatrixScalar,
00358          class DomainScalar,
00359          class RangeScalar>
00360 void
00361 upperTriSolveCscColMajor (
00362   const Ordinal numRows,
00363   const Ordinal numCols,
00364   const Ordinal numVecs,
00365   RangeScalar* const X,
00366   const Ordinal colStrideX,
00367   const  size_t* const ptr,
00368   const Ordinal* const ind,
00369   const MatrixScalar* const val,
00370   const DomainScalar* const Y,
00371   const Ordinal colStrideY)
00372 {
00373   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00374 
00375   for (Ordinal j = 0; j < numVecs; ++j) {
00376     for (Ordinal i = 0; i < numRows; ++i) {
00377       X[i + j*colStrideX] = Y[i + j*colStrideY];
00378     }
00379   }
00380 
00381   for (Ordinal c = numCols-1; c >= 0; --c) {
00382     MatrixScalar A_cc = STS::zero ();
00383     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00384       const Ordinal r = ind[k];
00385       MatrixScalar A_rc = val[k];
00386       if (r == c) {
00387         A_cc = A_cc + A_rc;
00388       } else {
00389         for (Ordinal j = 0; j < numVecs; ++j) {
00390           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00391         }
00392       }
00393       for (Ordinal j = 0; j < numVecs; ++j) {
00394         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
00395       }
00396     }
00397   }
00398 }
00399 
00400 template<class Ordinal,
00401          class MatrixScalar,
00402          class DomainScalar,
00403          class RangeScalar>
00404 void
00405 upperTriSolveCsrColMajor (
00406   const Ordinal numRows,
00407   const Ordinal numCols,
00408   const Ordinal numVecs,
00409   RangeScalar* const X,
00410   const Ordinal colStrideX,
00411   const  size_t* const ptr,
00412   const Ordinal* const ind,
00413   const MatrixScalar* const val,
00414   const DomainScalar* const Y,
00415   const Ordinal colStrideY)
00416 {
00417   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00418 
00419   for (Ordinal r = numRows-1; r >= 0; --r) {
00420     for (Ordinal j = 0; j < numVecs; ++j) {
00421       X[r + j*colStrideX] = Y[r + j*colStrideY];
00422     }
00423     // We assume the diagonal entry is first in the row.
00424     const MatrixScalar A_rr = val[ptr[r]];
00425     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
00426       const MatrixScalar A_rc = val[k];
00427       const Ordinal c = ind[k];
00428       for (Ordinal j = 0; j < numVecs; ++j) {
00429         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00430       }
00431     }
00432     for (Ordinal j = 0; j < numVecs; ++j) {
00433       X[r + j*colStrideX] = X[r + j*colStrideX] / A_rr;
00434     }
00435   }
00436 }
00437 
00438 template<class Ordinal,
00439          class MatrixScalar,
00440          class DomainScalar,
00441          class RangeScalar>
00442 void
00443 upperTriSolveCscRowMajor (
00444   const Ordinal numRows,
00445   const Ordinal numCols,
00446   const Ordinal numVecs,
00447   RangeScalar* const X,
00448   const Ordinal rowStrideX,
00449   const  size_t* const ptr,
00450   const Ordinal* const ind,
00451   const MatrixScalar* const val,
00452   const DomainScalar* const Y,
00453   const Ordinal rowStrideY)
00454 {
00455   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00456 
00457   for (Ordinal i = 0; i < numRows; ++i) {
00458     for (Ordinal j = 0; j < numVecs; ++j) {
00459       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
00460     }
00461   }
00462 
00463   for (Ordinal c = numCols-1; c >= 0; --c) {
00464     MatrixScalar A_cc = STS::zero ();
00465     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00466       const Ordinal r = ind[k];
00467       MatrixScalar A_rc = val[k];
00468       if (r == c) {
00469         A_cc = A_cc + A_rc;
00470       } else {
00471         for (Ordinal j = 0; j < numVecs; ++j) {
00472           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00473         }
00474       }
00475       for (Ordinal j = 0; j < numVecs; ++j) {
00476         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
00477       }
00478     }
00479   }
00480 }
00481 
00482 template<class Ordinal,
00483          class MatrixScalar,
00484          class DomainScalar,
00485          class RangeScalar>
00486 void
00487 upperTriSolveCsrRowMajor (
00488   const Ordinal numRows,
00489   const Ordinal numCols,
00490   const Ordinal numVecs,
00491   RangeScalar* const X,
00492   const Ordinal rowStrideX,
00493   const  size_t* const ptr,
00494   const Ordinal* const ind,
00495   const MatrixScalar* const val,
00496   const DomainScalar* const Y,
00497   const Ordinal rowStrideY)
00498 {
00499   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00500 
00501   for (Ordinal r = numRows-1; r >= 0; --r) {
00502     for (Ordinal j = 0; j < numVecs; ++j) {
00503       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
00504     }
00505     // We assume the diagonal entry is first in the row.
00506     const MatrixScalar A_rr = val[ptr[r]];
00507     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
00508       const MatrixScalar A_rc = val[k];
00509       const Ordinal c = ind[k];
00510       for (Ordinal j = 0; j < numVecs; ++j) {
00511         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00512       }
00513     }
00514     for (Ordinal j = 0; j < numVecs; ++j) {
00515       X[r*rowStrideX + j] = X[r*rowStrideX + j] / A_rr;
00516     }
00517   }
00518 }
00519 
00520 template<class Ordinal,
00521          class MatrixScalar,
00522          class DomainScalar,
00523          class RangeScalar>
00524 void
00525 upperTriSolveCscColMajorUnitDiag (
00526   const Ordinal numRows,
00527   const Ordinal numCols,
00528   const Ordinal numVecs,
00529   RangeScalar* const X,
00530   const Ordinal colStrideX,
00531   const  size_t* const ptr,
00532   const Ordinal* const ind,
00533   const MatrixScalar* const val,
00534   const DomainScalar* const Y,
00535   const Ordinal colStrideY)
00536 {
00537   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00538 
00539   for (Ordinal j = 0; j < numVecs; ++j) {
00540     for (Ordinal i = 0; i < numRows; ++i) {
00541       X[i + j*colStrideX] = Y[i + j*colStrideY];
00542     }
00543   }
00544 
00545   for (Ordinal c = numCols-1; c >= 0; --c) {
00546     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00547       const Ordinal r = ind[k];
00548       MatrixScalar A_rc = val[k];
00549       for (Ordinal j = 0; j < numVecs; ++j) {
00550         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00551       }
00552     }
00553   }
00554 }
00555 
00556 template<class Ordinal,
00557          class MatrixScalar,
00558          class DomainScalar,
00559          class RangeScalar>
00560 void
00561 upperTriSolveCsrColMajorUnitDiag (
00562   const Ordinal numRows,
00563   const Ordinal numCols,
00564   const Ordinal numVecs,
00565   RangeScalar* const X,
00566   const Ordinal colStrideX,
00567   const  size_t* const ptr,
00568   const Ordinal* const ind,
00569   const MatrixScalar* const val,
00570   const DomainScalar* const Y,
00571   const Ordinal colStrideY)
00572 {
00573   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00574 
00575   for (Ordinal r = numRows-1; r >= 0; --r) {
00576     for (Ordinal j = 0; j < numVecs; ++j) {
00577       X[r + j*colStrideX] = Y[r + j*colStrideY];
00578     }
00579     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
00580       const MatrixScalar A_rc = val[k];
00581       const Ordinal c = ind[k];
00582       for (Ordinal j = 0; j < numVecs; ++j) {
00583         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00584       }
00585     }
00586   }
00587 }
00588 
00589 template<class Ordinal,
00590          class MatrixScalar,
00591          class DomainScalar,
00592          class RangeScalar>
00593 void
00594 upperTriSolveCscRowMajorUnitDiag (
00595   const Ordinal numRows,
00596   const Ordinal numCols,
00597   const Ordinal numVecs,
00598   RangeScalar* const X,
00599   const Ordinal rowStrideX,
00600   const  size_t* const ptr,
00601   const Ordinal* const ind,
00602   const MatrixScalar* const val,
00603   const DomainScalar* const Y,
00604   const Ordinal rowStrideY)
00605 {
00606   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00607 
00608   for (Ordinal i = 0; i < numRows; ++i) {
00609     for (Ordinal j = 0; j < numVecs; ++j) {
00610       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
00611     }
00612   }
00613 
00614   for (Ordinal c = numCols-1; c >= 0; --c) {
00615     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00616       const Ordinal r = ind[k];
00617       MatrixScalar A_rc = val[k];
00618       for (Ordinal j = 0; j < numVecs; ++j) {
00619         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00620       }
00621     }
00622   }
00623 }
00624 
00625 template<class Ordinal,
00626          class MatrixScalar,
00627          class DomainScalar,
00628          class RangeScalar>
00629 void
00630 upperTriSolveCsrRowMajorUnitDiag (
00631   const Ordinal numRows,
00632   const Ordinal numCols,
00633   const Ordinal numVecs,
00634   RangeScalar* const X,
00635   const Ordinal rowStrideX,
00636   const  size_t* const ptr,
00637   const Ordinal* const ind,
00638   const MatrixScalar* const val,
00639   const DomainScalar* const Y,
00640   const Ordinal rowStrideY)
00641 {
00642   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00643 
00644   for (Ordinal r = numRows-1; r >= 0; --r) {
00645     for (Ordinal j = 0; j < numVecs; ++j) {
00646       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
00647     }
00648     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
00649       const MatrixScalar A_rc = val[k];
00650       const Ordinal c = ind[k];
00651       for (Ordinal j = 0; j < numVecs; ++j) {
00652         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00653       }
00654     }
00655   }
00656 }
00657 
00658 template<class Ordinal,
00659          class MatrixScalar,
00660          class RangeScalar>
00661 void
00662 lowerTriSolveCscColMajorInPlace (
00663   const Ordinal numRows,
00664   const Ordinal numCols,
00665   const Ordinal numVecs,
00666   RangeScalar* const X,
00667   const Ordinal colStrideX,
00668   const  size_t* const ptr,
00669   const Ordinal* const ind,
00670   const MatrixScalar* const val)
00671 {
00672   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00673 
00674   for (Ordinal c = 0; c < numCols; ++c) {
00675     MatrixScalar A_cc = STS::zero ();
00676     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00677       const Ordinal r = ind[k];
00678       MatrixScalar A_rc = val[k];
00679       if (r == c) {
00680         A_cc = A_cc + A_rc;
00681       } else {
00682         for (Ordinal j = 0; j < numVecs; ++j) {
00683           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00684         }
00685       }
00686       for (Ordinal j = 0; j < numVecs; ++j) {
00687         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
00688       }
00689     }
00690   }
00691 }
00692 
00693 template<class Ordinal,
00694          class MatrixScalar,
00695          class RangeScalar>
00696 void
00697 lowerTriSolveCscRowMajorInPlace (
00698   const Ordinal numRows,
00699   const Ordinal numCols,
00700   const Ordinal numVecs,
00701   RangeScalar* const X,
00702   const Ordinal rowStrideX,
00703   const  size_t* const ptr,
00704   const Ordinal* const ind,
00705   const MatrixScalar* const val)
00706 {
00707   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00708 
00709   for (Ordinal c = 0; c < numCols; ++c) {
00710     MatrixScalar A_cc = STS::zero ();
00711     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00712       const Ordinal r = ind[k];
00713       MatrixScalar A_rc = val[k];
00714       if (r == c) {
00715         A_cc = A_cc + A_rc;
00716       } else {
00717         for (Ordinal j = 0; j < numVecs; ++j) {
00718           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00719         }
00720       }
00721       for (Ordinal j = 0; j < numVecs; ++j) {
00722         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
00723       }
00724     }
00725   }
00726 }
00727 
00728 template<class Ordinal,
00729          class MatrixScalar,
00730          class RangeScalar>
00731 void
00732 lowerTriSolveCscColMajorUnitDiagInPlace (
00733   const Ordinal numRows,
00734   const Ordinal numCols,
00735   const Ordinal numVecs,
00736   RangeScalar* const X,
00737   const Ordinal colStrideX,
00738   const  size_t* const ptr,
00739   const Ordinal* const ind,
00740   const MatrixScalar* const val)
00741 {
00742   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00743 
00744   for (Ordinal c = 0; c < numCols; ++c) {
00745     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00746       const Ordinal r = ind[k];
00747       MatrixScalar A_rc = val[k];
00748       for (Ordinal j = 0; j < numVecs; ++j) {
00749         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00750       }
00751     }
00752   }
00753 }
00754 
00755 template<class Ordinal,
00756          class MatrixScalar,
00757          class RangeScalar>
00758 void
00759 lowerTriSolveCscRowMajorUnitDiagInPlace (
00760   const Ordinal numRows,
00761   const Ordinal numCols,
00762   const Ordinal numVecs,
00763   RangeScalar* const X,
00764   const Ordinal rowStrideX,
00765   const  size_t* const ptr,
00766   const Ordinal* const ind,
00767   const MatrixScalar* const val)
00768 {
00769   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00770 
00771   for (Ordinal c = 0; c < numCols; ++c) {
00772     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00773       const Ordinal r = ind[k];
00774       MatrixScalar A_rc = val[k];
00775       for (Ordinal j = 0; j < numVecs; ++j) {
00776         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00777       }
00778     }
00779   }
00780 }
00781 
00782 template<class Ordinal,
00783          class MatrixScalar,
00784          class RangeScalar>
00785 void
00786 upperTriSolveCscColMajorInPlace (
00787   const Ordinal numRows,
00788   const Ordinal numCols,
00789   const Ordinal numVecs,
00790   RangeScalar* const X,
00791   const Ordinal colStrideX,
00792   const  size_t* const ptr,
00793   const Ordinal* const ind,
00794   const MatrixScalar* const val)
00795 {
00796   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00797 
00798   for (Ordinal c = numCols-1; c >= 0; --c) {
00799     MatrixScalar A_cc = STS::zero ();
00800     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00801       const Ordinal r = ind[k];
00802       MatrixScalar A_rc = val[k];
00803       if (r == c) {
00804         A_cc = A_cc + A_rc;
00805       } else {
00806         for (Ordinal j = 0; j < numVecs; ++j) {
00807           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00808         }
00809       }
00810       for (Ordinal j = 0; j < numVecs; ++j) {
00811         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
00812       }
00813     }
00814   }
00815 }
00816 
00817 template<class Ordinal,
00818          class MatrixScalar,
00819          class RangeScalar>
00820 void
00821 upperTriSolveCscRowMajorInPlace (
00822   const Ordinal numRows,
00823   const Ordinal numCols,
00824   const Ordinal numVecs,
00825   RangeScalar* const X,
00826   const Ordinal rowStrideX,
00827   const  size_t* const ptr,
00828   const Ordinal* const ind,
00829   const MatrixScalar* const val)
00830 {
00831   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00832 
00833   for (Ordinal c = numCols-1; c >= 0; --c) {
00834     MatrixScalar A_cc = STS::zero ();
00835     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00836       const Ordinal r = ind[k];
00837       MatrixScalar A_rc = val[k];
00838       if (r == c) {
00839         A_cc = A_cc + A_rc;
00840       } else {
00841         for (Ordinal j = 0; j < numVecs; ++j) {
00842           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00843         }
00844       }
00845       for (Ordinal j = 0; j < numVecs; ++j) {
00846         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
00847       }
00848     }
00849   }
00850 }
00851 
00852 template<class Ordinal,
00853          class MatrixScalar,
00854          class RangeScalar>
00855 void
00856 upperTriSolveCscColMajorUnitDiagInPlace (
00857   const Ordinal numRows,
00858   const Ordinal numCols,
00859   const Ordinal numVecs,
00860   RangeScalar* const X,
00861   const Ordinal colStrideX,
00862   const  size_t* const ptr,
00863   const Ordinal* const ind,
00864   const MatrixScalar* const val)
00865 {
00866   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00867 
00868   for (Ordinal c = numCols-1; c >= 0; --c) {
00869     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00870       const Ordinal r = ind[k];
00871       MatrixScalar A_rc = val[k];
00872       for (Ordinal j = 0; j < numVecs; ++j) {
00873         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00874       }
00875     }
00876   }
00877 }
00878 
00879 template<class Ordinal,
00880          class MatrixScalar,
00881          class RangeScalar>
00882 void
00883 upperTriSolveCscRowMajorUnitDiagInPlace (
00884   const Ordinal numRows,
00885   const Ordinal numCols,
00886   const Ordinal numVecs,
00887   RangeScalar* const X,
00888   const Ordinal rowStrideX,
00889   const  size_t* const ptr,
00890   const Ordinal* const ind,
00891   const MatrixScalar* const val)
00892 {
00893   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00894 
00895   for (Ordinal c = numCols-1; c >= 0; --c) {
00896     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00897       const Ordinal r = ind[k];
00898       MatrixScalar A_rc = val[k];
00899       for (Ordinal j = 0; j < numVecs; ++j) {
00900         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
00901       }
00902     }
00903   }
00904 }
00905 
00906 template<class Ordinal,
00907          class MatrixScalar,
00908          class DomainScalar,
00909          class RangeScalar>
00910 void
00911 lowerTriSolveCscColMajorConj (
00912   const Ordinal numRows,
00913   const Ordinal numCols,
00914   const Ordinal numVecs,
00915   RangeScalar* const X,
00916   const Ordinal colStrideX,
00917   const  size_t* const ptr,
00918   const Ordinal* const ind,
00919   const MatrixScalar* const val,
00920   const DomainScalar* const Y,
00921   const Ordinal colStrideY)
00922 {
00923   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00924 
00925   for (Ordinal j = 0; j < numVecs; ++j) {
00926     for (Ordinal i = 0; i < numRows; ++i) {
00927       X[i + j*colStrideX] = Y[i + j*colStrideY];
00928     }
00929   }
00930 
00931   for (Ordinal c = 0; c < numCols; ++c) {
00932     MatrixScalar A_cc = STS::zero ();
00933     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
00934       const Ordinal r = ind[k];
00935       MatrixScalar A_rc = STS::conjugate (val[k]);
00936       if (r == c) {
00937         A_cc = A_cc + A_rc;
00938       } else {
00939         for (Ordinal j = 0; j < numVecs; ++j) {
00940           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00941         }
00942       }
00943       for (Ordinal j = 0; j < numVecs; ++j) {
00944         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
00945       }
00946     }
00947   }
00948 }
00949 
00950 template<class Ordinal,
00951          class MatrixScalar,
00952          class DomainScalar,
00953          class RangeScalar>
00954 void
00955 lowerTriSolveCsrColMajorConj (
00956   const Ordinal numRows,
00957   const Ordinal numCols,
00958   const Ordinal numVecs,
00959   RangeScalar* const X,
00960   const Ordinal colStrideX,
00961   const  size_t* const ptr,
00962   const Ordinal* const ind,
00963   const MatrixScalar* const val,
00964   const DomainScalar* const Y,
00965   const Ordinal colStrideY)
00966 {
00967   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
00968 
00969   for (Ordinal r = 0; r < numRows; ++r) {
00970     for (Ordinal j = 0; j < numVecs; ++j) {
00971       X[r + j*colStrideX] = Y[r + j*colStrideY];
00972     }
00973     // We assume the diagonal entry is first in the row.
00974     const MatrixScalar A_rr = STS::conjugate (val[ptr[r]]);
00975     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
00976       const MatrixScalar A_rc = STS::conjugate (val[k]);
00977       const Ordinal c = ind[k];
00978       for (Ordinal j = 0; j < numVecs; ++j) {
00979         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
00980       }
00981     }
00982     for (Ordinal j = 0; j < numVecs; ++j) {
00983       X[r + j*colStrideX] = X[r + j*colStrideX] / A_rr;
00984     }
00985   }
00986 }
00987 
00988 template<class Ordinal,
00989          class MatrixScalar,
00990          class DomainScalar,
00991          class RangeScalar>
00992 void
00993 lowerTriSolveCscRowMajorConj (
00994   const Ordinal numRows,
00995   const Ordinal numCols,
00996   const Ordinal numVecs,
00997   RangeScalar* const X,
00998   const Ordinal rowStrideX,
00999   const  size_t* const ptr,
01000   const Ordinal* const ind,
01001   const MatrixScalar* const val,
01002   const DomainScalar* const Y,
01003   const Ordinal rowStrideY)
01004 {
01005   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01006 
01007   for (Ordinal i = 0; i < numRows; ++i) {
01008     for (Ordinal j = 0; j < numVecs; ++j) {
01009       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
01010     }
01011   }
01012 
01013   for (Ordinal c = 0; c < numCols; ++c) {
01014     MatrixScalar A_cc = STS::zero ();
01015     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01016       const Ordinal r = ind[k];
01017       MatrixScalar A_rc = STS::conjugate (val[k]);
01018       if (r == c) {
01019         A_cc = A_cc + A_rc;
01020       } else {
01021         for (Ordinal j = 0; j < numVecs; ++j) {
01022           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01023         }
01024       }
01025       for (Ordinal j = 0; j < numVecs; ++j) {
01026         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
01027       }
01028     }
01029   }
01030 }
01031 
01032 template<class Ordinal,
01033          class MatrixScalar,
01034          class DomainScalar,
01035          class RangeScalar>
01036 void
01037 lowerTriSolveCsrRowMajorConj (
01038   const Ordinal numRows,
01039   const Ordinal numCols,
01040   const Ordinal numVecs,
01041   RangeScalar* const X,
01042   const Ordinal rowStrideX,
01043   const  size_t* const ptr,
01044   const Ordinal* const ind,
01045   const MatrixScalar* const val,
01046   const DomainScalar* const Y,
01047   const Ordinal rowStrideY)
01048 {
01049   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01050 
01051   for (Ordinal r = 0; r < numRows; ++r) {
01052     for (Ordinal j = 0; j < numVecs; ++j) {
01053       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
01054     }
01055     // We assume the diagonal entry is first in the row.
01056     const MatrixScalar A_rr = STS::conjugate (val[ptr[r]]);
01057     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
01058       const MatrixScalar A_rc = STS::conjugate (val[k]);
01059       const Ordinal c = ind[k];
01060       for (Ordinal j = 0; j < numVecs; ++j) {
01061         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01062       }
01063     }
01064     for (Ordinal j = 0; j < numVecs; ++j) {
01065       X[r*rowStrideX + j] = X[r*rowStrideX + j] / A_rr;
01066     }
01067   }
01068 }
01069 
01070 template<class Ordinal,
01071          class MatrixScalar,
01072          class DomainScalar,
01073          class RangeScalar>
01074 void
01075 lowerTriSolveCscColMajorUnitDiagConj (
01076   const Ordinal numRows,
01077   const Ordinal numCols,
01078   const Ordinal numVecs,
01079   RangeScalar* const X,
01080   const Ordinal colStrideX,
01081   const  size_t* const ptr,
01082   const Ordinal* const ind,
01083   const MatrixScalar* const val,
01084   const DomainScalar* const Y,
01085   const Ordinal colStrideY)
01086 {
01087   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01088 
01089   for (Ordinal j = 0; j < numVecs; ++j) {
01090     for (Ordinal i = 0; i < numRows; ++i) {
01091       X[i + j*colStrideX] = Y[i + j*colStrideY];
01092     }
01093   }
01094 
01095   for (Ordinal c = 0; c < numCols; ++c) {
01096     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01097       const Ordinal r = ind[k];
01098       MatrixScalar A_rc = STS::conjugate (val[k]);
01099       for (Ordinal j = 0; j < numVecs; ++j) {
01100         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01101       }
01102     }
01103   }
01104 }
01105 
01106 template<class Ordinal,
01107          class MatrixScalar,
01108          class DomainScalar,
01109          class RangeScalar>
01110 void
01111 lowerTriSolveCsrColMajorUnitDiagConj (
01112   const Ordinal numRows,
01113   const Ordinal numCols,
01114   const Ordinal numVecs,
01115   RangeScalar* const X,
01116   const Ordinal colStrideX,
01117   const  size_t* const ptr,
01118   const Ordinal* const ind,
01119   const MatrixScalar* const val,
01120   const DomainScalar* const Y,
01121   const Ordinal colStrideY)
01122 {
01123   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01124 
01125   for (Ordinal r = 0; r < numRows; ++r) {
01126     for (Ordinal j = 0; j < numVecs; ++j) {
01127       X[r + j*colStrideX] = Y[r + j*colStrideY];
01128     }
01129     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
01130       const MatrixScalar A_rc = STS::conjugate (val[k]);
01131       const Ordinal c = ind[k];
01132       for (Ordinal j = 0; j < numVecs; ++j) {
01133         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01134       }
01135     }
01136   }
01137 }
01138 
01139 template<class Ordinal,
01140          class MatrixScalar,
01141          class DomainScalar,
01142          class RangeScalar>
01143 void
01144 lowerTriSolveCscRowMajorUnitDiagConj (
01145   const Ordinal numRows,
01146   const Ordinal numCols,
01147   const Ordinal numVecs,
01148   RangeScalar* const X,
01149   const Ordinal rowStrideX,
01150   const  size_t* const ptr,
01151   const Ordinal* const ind,
01152   const MatrixScalar* const val,
01153   const DomainScalar* const Y,
01154   const Ordinal rowStrideY)
01155 {
01156   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01157 
01158   for (Ordinal i = 0; i < numRows; ++i) {
01159     for (Ordinal j = 0; j < numVecs; ++j) {
01160       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
01161     }
01162   }
01163 
01164   for (Ordinal c = 0; c < numCols; ++c) {
01165     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01166       const Ordinal r = ind[k];
01167       MatrixScalar A_rc = STS::conjugate (val[k]);
01168       for (Ordinal j = 0; j < numVecs; ++j) {
01169         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01170       }
01171     }
01172   }
01173 }
01174 
01175 template<class Ordinal,
01176          class MatrixScalar,
01177          class DomainScalar,
01178          class RangeScalar>
01179 void
01180 lowerTriSolveCsrRowMajorUnitDiagConj (
01181   const Ordinal numRows,
01182   const Ordinal numCols,
01183   const Ordinal numVecs,
01184   RangeScalar* const X,
01185   const Ordinal rowStrideX,
01186   const  size_t* const ptr,
01187   const Ordinal* const ind,
01188   const MatrixScalar* const val,
01189   const DomainScalar* const Y,
01190   const Ordinal rowStrideY)
01191 {
01192   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01193 
01194   for (Ordinal r = 0; r < numRows; ++r) {
01195     for (Ordinal j = 0; j < numVecs; ++j) {
01196       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
01197     }
01198     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
01199       const MatrixScalar A_rc = STS::conjugate (val[k]);
01200       const Ordinal c = ind[k];
01201       for (Ordinal j = 0; j < numVecs; ++j) {
01202         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01203       }
01204     }
01205   }
01206 }
01207 
01208 template<class Ordinal,
01209          class MatrixScalar,
01210          class DomainScalar,
01211          class RangeScalar>
01212 void
01213 upperTriSolveCscColMajorConj (
01214   const Ordinal numRows,
01215   const Ordinal numCols,
01216   const Ordinal numVecs,
01217   RangeScalar* const X,
01218   const Ordinal colStrideX,
01219   const  size_t* const ptr,
01220   const Ordinal* const ind,
01221   const MatrixScalar* const val,
01222   const DomainScalar* const Y,
01223   const Ordinal colStrideY)
01224 {
01225   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01226 
01227   for (Ordinal j = 0; j < numVecs; ++j) {
01228     for (Ordinal i = 0; i < numRows; ++i) {
01229       X[i + j*colStrideX] = Y[i + j*colStrideY];
01230     }
01231   }
01232 
01233   for (Ordinal c = numCols-1; c >= 0; --c) {
01234     MatrixScalar A_cc = STS::zero ();
01235     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01236       const Ordinal r = ind[k];
01237       MatrixScalar A_rc = STS::conjugate (val[k]);
01238       if (r == c) {
01239         A_cc = A_cc + A_rc;
01240       } else {
01241         for (Ordinal j = 0; j < numVecs; ++j) {
01242           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01243         }
01244       }
01245       for (Ordinal j = 0; j < numVecs; ++j) {
01246         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
01247       }
01248     }
01249   }
01250 }
01251 
01252 template<class Ordinal,
01253          class MatrixScalar,
01254          class DomainScalar,
01255          class RangeScalar>
01256 void
01257 upperTriSolveCsrColMajorConj (
01258   const Ordinal numRows,
01259   const Ordinal numCols,
01260   const Ordinal numVecs,
01261   RangeScalar* const X,
01262   const Ordinal colStrideX,
01263   const  size_t* const ptr,
01264   const Ordinal* const ind,
01265   const MatrixScalar* const val,
01266   const DomainScalar* const Y,
01267   const Ordinal colStrideY)
01268 {
01269   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01270 
01271   for (Ordinal r = numRows-1; r >= 0; --r) {
01272     for (Ordinal j = 0; j < numVecs; ++j) {
01273       X[r + j*colStrideX] = Y[r + j*colStrideY];
01274     }
01275     // We assume the diagonal entry is first in the row.
01276     const MatrixScalar A_rr = STS::conjugate (val[ptr[r]]);
01277     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
01278       const MatrixScalar A_rc = STS::conjugate (val[k]);
01279       const Ordinal c = ind[k];
01280       for (Ordinal j = 0; j < numVecs; ++j) {
01281         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01282       }
01283     }
01284     for (Ordinal j = 0; j < numVecs; ++j) {
01285       X[r + j*colStrideX] = X[r + j*colStrideX] / A_rr;
01286     }
01287   }
01288 }
01289 
01290 template<class Ordinal,
01291          class MatrixScalar,
01292          class DomainScalar,
01293          class RangeScalar>
01294 void
01295 upperTriSolveCscRowMajorConj (
01296   const Ordinal numRows,
01297   const Ordinal numCols,
01298   const Ordinal numVecs,
01299   RangeScalar* const X,
01300   const Ordinal rowStrideX,
01301   const  size_t* const ptr,
01302   const Ordinal* const ind,
01303   const MatrixScalar* const val,
01304   const DomainScalar* const Y,
01305   const Ordinal rowStrideY)
01306 {
01307   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01308 
01309   for (Ordinal i = 0; i < numRows; ++i) {
01310     for (Ordinal j = 0; j < numVecs; ++j) {
01311       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
01312     }
01313   }
01314 
01315   for (Ordinal c = numCols-1; c >= 0; --c) {
01316     MatrixScalar A_cc = STS::zero ();
01317     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01318       const Ordinal r = ind[k];
01319       MatrixScalar A_rc = STS::conjugate (val[k]);
01320       if (r == c) {
01321         A_cc = A_cc + A_rc;
01322       } else {
01323         for (Ordinal j = 0; j < numVecs; ++j) {
01324           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01325         }
01326       }
01327       for (Ordinal j = 0; j < numVecs; ++j) {
01328         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
01329       }
01330     }
01331   }
01332 }
01333 
01334 template<class Ordinal,
01335          class MatrixScalar,
01336          class DomainScalar,
01337          class RangeScalar>
01338 void
01339 upperTriSolveCsrRowMajorConj (
01340   const Ordinal numRows,
01341   const Ordinal numCols,
01342   const Ordinal numVecs,
01343   RangeScalar* const X,
01344   const Ordinal rowStrideX,
01345   const  size_t* const ptr,
01346   const Ordinal* const ind,
01347   const MatrixScalar* const val,
01348   const DomainScalar* const Y,
01349   const Ordinal rowStrideY)
01350 {
01351   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01352 
01353   for (Ordinal r = numRows-1; r >= 0; --r) {
01354     for (Ordinal j = 0; j < numVecs; ++j) {
01355       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
01356     }
01357     // We assume the diagonal entry is first in the row.
01358     const MatrixScalar A_rr = STS::conjugate (val[ptr[r]]);
01359     for (size_t k = ptr[r]+1; k < ptr[r+1]; ++k) {
01360       const MatrixScalar A_rc = STS::conjugate (val[k]);
01361       const Ordinal c = ind[k];
01362       for (Ordinal j = 0; j < numVecs; ++j) {
01363         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01364       }
01365     }
01366     for (Ordinal j = 0; j < numVecs; ++j) {
01367       X[r*rowStrideX + j] = X[r*rowStrideX + j] / A_rr;
01368     }
01369   }
01370 }
01371 
01372 template<class Ordinal,
01373          class MatrixScalar,
01374          class DomainScalar,
01375          class RangeScalar>
01376 void
01377 upperTriSolveCscColMajorUnitDiagConj (
01378   const Ordinal numRows,
01379   const Ordinal numCols,
01380   const Ordinal numVecs,
01381   RangeScalar* const X,
01382   const Ordinal colStrideX,
01383   const  size_t* const ptr,
01384   const Ordinal* const ind,
01385   const MatrixScalar* const val,
01386   const DomainScalar* const Y,
01387   const Ordinal colStrideY)
01388 {
01389   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01390 
01391   for (Ordinal j = 0; j < numVecs; ++j) {
01392     for (Ordinal i = 0; i < numRows; ++i) {
01393       X[i + j*colStrideX] = Y[i + j*colStrideY];
01394     }
01395   }
01396 
01397   for (Ordinal c = numCols-1; c >= 0; --c) {
01398     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01399       const Ordinal r = ind[k];
01400       MatrixScalar A_rc = STS::conjugate (val[k]);
01401       for (Ordinal j = 0; j < numVecs; ++j) {
01402         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01403       }
01404     }
01405   }
01406 }
01407 
01408 template<class Ordinal,
01409          class MatrixScalar,
01410          class DomainScalar,
01411          class RangeScalar>
01412 void
01413 upperTriSolveCsrColMajorUnitDiagConj (
01414   const Ordinal numRows,
01415   const Ordinal numCols,
01416   const Ordinal numVecs,
01417   RangeScalar* const X,
01418   const Ordinal colStrideX,
01419   const  size_t* const ptr,
01420   const Ordinal* const ind,
01421   const MatrixScalar* const val,
01422   const DomainScalar* const Y,
01423   const Ordinal colStrideY)
01424 {
01425   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01426 
01427   for (Ordinal r = numRows-1; r >= 0; --r) {
01428     for (Ordinal j = 0; j < numVecs; ++j) {
01429       X[r + j*colStrideX] = Y[r + j*colStrideY];
01430     }
01431     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
01432       const MatrixScalar A_rc = STS::conjugate (val[k]);
01433       const Ordinal c = ind[k];
01434       for (Ordinal j = 0; j < numVecs; ++j) {
01435         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01436       }
01437     }
01438   }
01439 }
01440 
01441 template<class Ordinal,
01442          class MatrixScalar,
01443          class DomainScalar,
01444          class RangeScalar>
01445 void
01446 upperTriSolveCscRowMajorUnitDiagConj (
01447   const Ordinal numRows,
01448   const Ordinal numCols,
01449   const Ordinal numVecs,
01450   RangeScalar* const X,
01451   const Ordinal rowStrideX,
01452   const  size_t* const ptr,
01453   const Ordinal* const ind,
01454   const MatrixScalar* const val,
01455   const DomainScalar* const Y,
01456   const Ordinal rowStrideY)
01457 {
01458   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01459 
01460   for (Ordinal i = 0; i < numRows; ++i) {
01461     for (Ordinal j = 0; j < numVecs; ++j) {
01462       X[i*rowStrideX + j] = Y[i*rowStrideY + j];
01463     }
01464   }
01465 
01466   for (Ordinal c = numCols-1; c >= 0; --c) {
01467     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01468       const Ordinal r = ind[k];
01469       MatrixScalar A_rc = STS::conjugate (val[k]);
01470       for (Ordinal j = 0; j < numVecs; ++j) {
01471         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01472       }
01473     }
01474   }
01475 }
01476 
01477 template<class Ordinal,
01478          class MatrixScalar,
01479          class DomainScalar,
01480          class RangeScalar>
01481 void
01482 upperTriSolveCsrRowMajorUnitDiagConj (
01483   const Ordinal numRows,
01484   const Ordinal numCols,
01485   const Ordinal numVecs,
01486   RangeScalar* const X,
01487   const Ordinal rowStrideX,
01488   const  size_t* const ptr,
01489   const Ordinal* const ind,
01490   const MatrixScalar* const val,
01491   const DomainScalar* const Y,
01492   const Ordinal rowStrideY)
01493 {
01494   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01495 
01496   for (Ordinal r = numRows-1; r >= 0; --r) {
01497     for (Ordinal j = 0; j < numVecs; ++j) {
01498       X[r*rowStrideX + j] = Y[r*rowStrideY + j];
01499     }
01500     for (size_t k = ptr[r]; k < ptr[r+1]; ++k) {
01501       const MatrixScalar A_rc = STS::conjugate (val[k]);
01502       const Ordinal c = ind[k];
01503       for (Ordinal j = 0; j < numVecs; ++j) {
01504         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01505       }
01506     }
01507   }
01508 }
01509 
01510 template<class Ordinal,
01511          class MatrixScalar,
01512          class RangeScalar>
01513 void
01514 lowerTriSolveCscColMajorInPlaceConj (
01515   const Ordinal numRows,
01516   const Ordinal numCols,
01517   const Ordinal numVecs,
01518   RangeScalar* const X,
01519   const Ordinal colStrideX,
01520   const  size_t* const ptr,
01521   const Ordinal* const ind,
01522   const MatrixScalar* const val)
01523 {
01524   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01525 
01526   for (Ordinal c = 0; c < numCols; ++c) {
01527     MatrixScalar A_cc = STS::zero ();
01528     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01529       const Ordinal r = ind[k];
01530       MatrixScalar A_rc = STS::conjugate (val[k]);
01531       if (r == c) {
01532         A_cc = A_cc + A_rc;
01533       } else {
01534         for (Ordinal j = 0; j < numVecs; ++j) {
01535           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01536         }
01537       }
01538       for (Ordinal j = 0; j < numVecs; ++j) {
01539         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
01540       }
01541     }
01542   }
01543 }
01544 
01545 template<class Ordinal,
01546          class MatrixScalar,
01547          class RangeScalar>
01548 void
01549 lowerTriSolveCscRowMajorInPlaceConj (
01550   const Ordinal numRows,
01551   const Ordinal numCols,
01552   const Ordinal numVecs,
01553   RangeScalar* const X,
01554   const Ordinal rowStrideX,
01555   const  size_t* const ptr,
01556   const Ordinal* const ind,
01557   const MatrixScalar* const val)
01558 {
01559   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01560 
01561   for (Ordinal c = 0; c < numCols; ++c) {
01562     MatrixScalar A_cc = STS::zero ();
01563     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01564       const Ordinal r = ind[k];
01565       MatrixScalar A_rc = STS::conjugate (val[k]);
01566       if (r == c) {
01567         A_cc = A_cc + A_rc;
01568       } else {
01569         for (Ordinal j = 0; j < numVecs; ++j) {
01570           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01571         }
01572       }
01573       for (Ordinal j = 0; j < numVecs; ++j) {
01574         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
01575       }
01576     }
01577   }
01578 }
01579 
01580 template<class Ordinal,
01581          class MatrixScalar,
01582          class RangeScalar>
01583 void
01584 lowerTriSolveCscColMajorUnitDiagInPlaceConj (
01585   const Ordinal numRows,
01586   const Ordinal numCols,
01587   const Ordinal numVecs,
01588   RangeScalar* const X,
01589   const Ordinal colStrideX,
01590   const  size_t* const ptr,
01591   const Ordinal* const ind,
01592   const MatrixScalar* const val)
01593 {
01594   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01595 
01596   for (Ordinal c = 0; c < numCols; ++c) {
01597     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01598       const Ordinal r = ind[k];
01599       MatrixScalar A_rc = STS::conjugate (val[k]);
01600       for (Ordinal j = 0; j < numVecs; ++j) {
01601         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01602       }
01603     }
01604   }
01605 }
01606 
01607 template<class Ordinal,
01608          class MatrixScalar,
01609          class RangeScalar>
01610 void
01611 lowerTriSolveCscRowMajorUnitDiagInPlaceConj (
01612   const Ordinal numRows,
01613   const Ordinal numCols,
01614   const Ordinal numVecs,
01615   RangeScalar* const X,
01616   const Ordinal rowStrideX,
01617   const  size_t* const ptr,
01618   const Ordinal* const ind,
01619   const MatrixScalar* const val)
01620 {
01621   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01622 
01623   for (Ordinal c = 0; c < numCols; ++c) {
01624     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01625       const Ordinal r = ind[k];
01626       MatrixScalar A_rc = STS::conjugate (val[k]);
01627       for (Ordinal j = 0; j < numVecs; ++j) {
01628         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01629       }
01630     }
01631   }
01632 }
01633 
01634 template<class Ordinal,
01635          class MatrixScalar,
01636          class RangeScalar>
01637 void
01638 upperTriSolveCscColMajorInPlaceConj (
01639   const Ordinal numRows,
01640   const Ordinal numCols,
01641   const Ordinal numVecs,
01642   RangeScalar* const X,
01643   const Ordinal colStrideX,
01644   const  size_t* const ptr,
01645   const Ordinal* const ind,
01646   const MatrixScalar* const val)
01647 {
01648   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01649 
01650   for (Ordinal c = numCols-1; c >= 0; --c) {
01651     MatrixScalar A_cc = STS::zero ();
01652     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01653       const Ordinal r = ind[k];
01654       MatrixScalar A_rc = STS::conjugate (val[k]);
01655       if (r == c) {
01656         A_cc = A_cc + A_rc;
01657       } else {
01658         for (Ordinal j = 0; j < numVecs; ++j) {
01659           X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01660         }
01661       }
01662       for (Ordinal j = 0; j < numVecs; ++j) {
01663         X[c + j*colStrideX] = X[c + j*colStrideX] / A_cc;
01664       }
01665     }
01666   }
01667 }
01668 
01669 template<class Ordinal,
01670          class MatrixScalar,
01671          class RangeScalar>
01672 void
01673 upperTriSolveCscRowMajorInPlaceConj (
01674   const Ordinal numRows,
01675   const Ordinal numCols,
01676   const Ordinal numVecs,
01677   RangeScalar* const X,
01678   const Ordinal rowStrideX,
01679   const  size_t* const ptr,
01680   const Ordinal* const ind,
01681   const MatrixScalar* const val)
01682 {
01683   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01684 
01685   for (Ordinal c = numCols-1; c >= 0; --c) {
01686     MatrixScalar A_cc = STS::zero ();
01687     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01688       const Ordinal r = ind[k];
01689       MatrixScalar A_rc = STS::conjugate (val[k]);
01690       if (r == c) {
01691         A_cc = A_cc + A_rc;
01692       } else {
01693         for (Ordinal j = 0; j < numVecs; ++j) {
01694           X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01695         }
01696       }
01697       for (Ordinal j = 0; j < numVecs; ++j) {
01698         X[c*rowStrideX + j] = X[c*rowStrideX + j] / A_cc;
01699       }
01700     }
01701   }
01702 }
01703 
01704 template<class Ordinal,
01705          class MatrixScalar,
01706          class RangeScalar>
01707 void
01708 upperTriSolveCscColMajorUnitDiagInPlaceConj (
01709   const Ordinal numRows,
01710   const Ordinal numCols,
01711   const Ordinal numVecs,
01712   RangeScalar* const X,
01713   const Ordinal colStrideX,
01714   const  size_t* const ptr,
01715   const Ordinal* const ind,
01716   const MatrixScalar* const val)
01717 {
01718   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01719 
01720   for (Ordinal c = numCols-1; c >= 0; --c) {
01721     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01722       const Ordinal r = ind[k];
01723       MatrixScalar A_rc = STS::conjugate (val[k]);
01724       for (Ordinal j = 0; j < numVecs; ++j) {
01725         X[r + j*colStrideX] -= A_rc * X[c + j*colStrideX];
01726       }
01727     }
01728   }
01729 }
01730 
01731 template<class Ordinal,
01732          class MatrixScalar,
01733          class RangeScalar>
01734 void
01735 upperTriSolveCscRowMajorUnitDiagInPlaceConj (
01736   const Ordinal numRows,
01737   const Ordinal numCols,
01738   const Ordinal numVecs,
01739   RangeScalar* const X,
01740   const Ordinal rowStrideX,
01741   const  size_t* const ptr,
01742   const Ordinal* const ind,
01743   const MatrixScalar* const val)
01744 {
01745   typedef Teuchos::ScalarTraits<MatrixScalar> STS;
01746 
01747   for (Ordinal c = numCols-1; c >= 0; --c) {
01748     for (size_t k = ptr[c]; k < ptr[c+1]; ++k) {
01749       const Ordinal r = ind[k];
01750       MatrixScalar A_rc = STS::conjugate (val[k]);
01751       for (Ordinal j = 0; j < numVecs; ++j) {
01752         X[r*rowStrideX + j] -= A_rc * X[c*rowStrideX + j];
01753       }
01754     }
01755   }
01756 }
01757 
01758 } // namespace Raw
01759 } // namespace Kokkos
01760 
01761 #endif // #ifndef __Kokkos_Raw_SparseTriangularSolve_def_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends