Anasazi Version of the Day
Tsqr_ZLapack.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <Tsqr_Lapack.hpp>
00030 #include <Tsqr_Config.hpp>
00031 #include <complex>
00032 
00035 
00036 extern "C" void F77_BLAS_MANGLE(zlarnv, ZLARNV)
00037   (const int* const IDIST,
00038    int ISEED[],
00039    const int* const N,
00040    std::complex<double> X[]);
00041 
00042 extern "C" void F77_BLAS_MANGLE(zpotri, ZPOTRI)
00043   (const char* const UPLO,
00044    const int* const N,
00045    std::complex<double> A[],
00046    const int* const LDA,
00047    int* const INFO);
00048 
00049 extern "C" void F77_BLAS_MANGLE(zpotrf, ZPOTRF)
00050   (const char* const UPLO,
00051    const int* const N,
00052    std::complex<double> A[],
00053    const int* const LDA,
00054    int* const INFO);
00055 
00056 extern "C" void F77_BLAS_MANGLE(zpotrs, ZPOTRS)
00057   (const char* const UPLO,
00058    const int* const N,
00059    const int* const NRHS,
00060    const std::complex<double> A[],
00061    const int* const LDA,
00062    std::complex<double> B[],
00063    const int* const LDB,
00064    int* const INFO);
00065 
00066 #ifdef HAVE_LAPACK_ZLARFGP
00067 extern "C" void F77_BLAS_MANGLE(zlarfgp,ZLARFGP)
00068   (const int* const N,    // IN
00069    std::complex<double>* const ALPHA,   // IN/OUT
00070    std::complex<double> X[],            // IN/OUT
00071    const int* const INCX, // IN
00072    std::complex<double>* const TAU);    // OUT
00073 #else
00074 #  ifdef HAVE_LAPACK_ZLARFP
00075 extern "C" void F77_BLAS_MANGLE(zlarfp,ZLARFP)
00076   (const int* const N,    // IN
00077    std::complex<double>* const ALPHA,   // IN/OUT
00078    std::complex<double> X[],            // IN/OUT
00079    const int* const INCX, // IN
00080    std::complex<double>* const TAU);    // OUT
00081 #  else
00082 extern "C" void F77_BLAS_MANGLE(zlarfg,ZLARFG)
00083   (const int* const N,    // IN
00084    std::complex<double>* const ALPHA,   // IN/OUT
00085    std::complex<double> X[],            // IN/OUT
00086    const int* const INCX, // IN
00087    std::complex<double>* const TAU);    // OUT
00088 #  endif // HAVE_LAPACK_ZLARFP
00089 #endif // HAVE_LAPACK_ZLARFGP
00090 
00091 extern "C" void F77_BLAS_MANGLE(zgeqrf, ZGEQRF)
00092   (const int* const M,
00093    const int* const N,
00094    std::complex<double> A[],
00095    const int* const LDA,
00096    std::complex<double> TAU[],
00097    std::complex<double> WORK[],
00098    const int* const LWORK,
00099    int* const INFO);
00100 
00101 #ifdef HAVE_LAPACK_ZGEQRFP
00102 extern "C" void F77_BLAS_MANGLE(zgeqrfp, ZGEQRFP)
00103   (const int* const M,
00104    const int* const N,
00105    std::complex<double> A[],
00106    const int* const LDA,
00107    std::complex<double> TAU[],
00108    std::complex<double> WORK[],
00109    const int* const LWORK,
00110    int* const INFO);
00111 #endif // HAVE_LAPACK_ZGEQRFP
00112 
00113 extern "C" void F77_BLAS_MANGLE(zgeqr2, ZGEQR2)
00114   (const int* const M,
00115    const int* const N,
00116    std::complex<double> A[],
00117    const int* const LDA,
00118    std::complex<double> TAU[],
00119    std::complex<double> WORK[],
00120    int* const INFO);
00121 
00122 #ifdef HAVE_LAPACK_ZGEQR2P
00123 extern "C" void F77_BLAS_MANGLE(zgeqr2p, ZGEQR2P)
00124   (const int* const M,
00125    const int* const N,
00126    std::complex<double> A[],
00127    const int* const LDA,
00128    std::complex<double> TAU[],
00129    std::complex<double> WORK[],
00130    int* const INFO);
00131 #endif // HAVE_LAPACK_ZGEQR2P
00132 
00133 // In the complex case, Q is called UNitary rather than ORthogonal.
00134 // This is why we have ZUNGQR and CUNGQR, rather than ZORGQR and
00135 // CORGQR.  The interface is exactly the same as in the real case,
00136 // though, so our LAPACK::ORMQR(), etc. wrappers have the same name
00137 // for both the real and the complex cases.
00138 
00139 extern "C" void F77_BLAS_MANGLE(zungqr, ZUNGQR)
00140   (const int* const M,
00141    const int* const N,
00142    const int* const K,
00143    std::complex<double> A[],
00144    const int* const LDA,
00145    std::complex<double> TAU[],
00146    std::complex<double> WORK[],
00147    const int* const LWORK,
00148    int* const INFO);
00149 
00150 extern "C" void F77_BLAS_MANGLE(zunmqr, ZUNMQR)
00151   (const char* const SIDE,
00152    const char* const TRANS,
00153    const int* const M,
00154    const int* const N,
00155    const int* const K,
00156    const std::complex<double> A[],
00157    const int* const LDA,
00158    const std::complex<double> TAU[],
00159    std::complex<double> C[],
00160    const int* const LDC,
00161    std::complex<double> WORK[],
00162    const int* const LWORK,
00163    int* const INFO);
00164 
00165 extern "C" void F77_BLAS_MANGLE(zunm2r, ZUNM2R)
00166   (const char* const SIDE,
00167    const char* const TRANS,
00168    const int* const M,
00169    const int* const N,
00170    const int* const K,
00171    const std::complex<double> A[],
00172    const int* const LDA,
00173    const std::complex<double> TAU[],
00174    std::complex<double> C[],
00175    const int* const LDC,
00176    std::complex<double> WORK[],
00177    int* const INFO);
00178 
00179 extern "C" void F77_BLAS_MANGLE(zgesvd, ZGESVD) 
00180   (const char* const JOBU, 
00181    const char* const JOBVT, 
00182    const int* const M, 
00183    const int* const N, 
00184    std::complex<double> A[], 
00185    const int* const LDA,
00186    double S[], 
00187    std::complex<double> U[], 
00188    const int* const LDU, 
00189    std::complex<double> VT[], 
00190    const int* const LDVT, 
00191    std::complex<double> work[],
00192    const int* const LWORK,
00193    double RWORK[],
00194    int* const INFO);
00195 
00198 
00199 namespace TSQR {
00200 
00201   // If _GEQRFP is available, LAPACK::GEQRF() calls it.  If _LARFP is
00202   // available, LAPACK::GEQRF() calls _GEQRF, which uses _LARFP.
00203 #ifdef HAVE_LAPACK_ZGEQRFP
00204   template <>
00205   bool LAPACK<int, std::complex<double> >::QR_produces_R_factor_with_nonnegative_diagonal() { return true; }
00206 #else
00207 #  ifdef HAVE_LAPACK_ZLARFP
00208   template <>
00209   bool LAPACK<int, std::complex<double> >::QR_produces_R_factor_with_nonnegative_diagonal() { return true; }
00210 #  else
00211   template <>
00212   bool LAPACK<int, std::complex<double> >::QR_produces_R_factor_with_nonnegative_diagonal() { return false; }
00213 #  endif
00214 #endif
00215 
00217   // LARFP (implemented with _LARFGP if available, else with _LARFP if
00218   // available, else fall back to _LARFG)
00220   template <>
00221   void 
00222   LAPACK<int, std::complex<double> >::
00223   LARFP (const int n, 
00224    std::complex<double>& alpha, 
00225    std::complex<double> x[], 
00226    const int incx, 
00227    std::complex<double>& tau)
00228   {
00229 #ifdef HAVE_LAPACK_ZLARFGP
00230     F77_BLAS_MANGLE(zlarfgp,ZLARFGP) (&n, &alpha, x, &incx, &tau);
00231 #else // Don't HAVE_LAPACK_CLARFGP
00232 #  ifdef HAVE_LAPACK_ZLARFP
00233     F77_BLAS_MANGLE(zlarfp,ZLARFP) (&n, &alpha, x, &incx, &tau);
00234 #  else
00235     F77_BLAS_MANGLE(zlarfg,ZLARFG) (&n, &alpha, x, &incx, &tau);
00236 #  endif // HAVE_LAPACK_ZLARFP
00237 #endif // HAVE_LAPACK_ZLARFGP
00238   }
00239 
00241   // GEQRF (implemented with _GEQRFP if available, else fall back to _GEQRF)
00243   template <>
00244   void
00245   LAPACK<int, std::complex<double> >::
00246   GEQRF (const int m,
00247    const int n, 
00248    std::complex<double> A[],
00249    const int lda, 
00250    std::complex<double> tau[],
00251    std::complex<double> work[],
00252    const int lwork,
00253    int* const INFO)
00254   {
00255 #ifdef HAVE_LAPACK_ZGEQRFP
00256     F77_BLAS_MANGLE(zgeqrfp, ZGEQRFP) (&m, &n, A, &lda, tau, work, &lwork, INFO);
00257 #else
00258     F77_BLAS_MANGLE(zgeqrf, ZGEQRF) (&m, &n, A, &lda, tau, work, &lwork, INFO);
00259 #endif // HAVE_LAPACK_ZGEQRFP
00260   }
00261 
00263   // GEQR2 (implemented with _GEQR2P if available, else fall back to _GEQR2)
00265   template <>
00266   void
00267   LAPACK<int, std::complex<double> >::
00268   GEQR2 (const int m,
00269    const int n, 
00270    std::complex<double> A[],
00271    const int lda, 
00272    std::complex<double> tau[],
00273    std::complex<double> work[],
00274    int* const INFO)
00275   {
00276 #ifdef HAVE_LAPACK_ZGEQR2P
00277     F77_BLAS_MANGLE(zgeqr2p, ZGEQR2P) (&m, &n, A, &lda, tau, work, INFO);
00278 #else
00279     F77_BLAS_MANGLE(zgeqr2, ZGEQR2) (&m, &n, A, &lda, tau, work, INFO);
00280 #endif // HAVE_LAPACK_ZGEQR2P
00281   }
00282 
00283   template <>
00284   void
00285   LAPACK<int, std::complex<double> >::
00286   ORMQR (const char* const side,
00287    const char* const trans,
00288    const int m,
00289    const int n,
00290    const int k,
00291    const std::complex<double> A[],
00292    const int lda,
00293    const std::complex<double> tau[],
00294    std::complex<double> C[],
00295    const int ldc,
00296    std::complex<double> work[],
00297    const int lwork,
00298    int* const INFO)
00299   {
00300     F77_BLAS_MANGLE(zunmqr, ZUNMQR) 
00301       (side, trans, &m, &n, &k, A, &lda, tau, C, &ldc, work, &lwork, INFO);
00302   }
00303 
00304   template <>
00305   void
00306   LAPACK<int, std::complex<double> >::
00307   ORM2R (const char* const side,
00308    const char* const trans,
00309    const int m,
00310    const int n,
00311    const int k,
00312    const std::complex<double> A[],
00313    const int lda,
00314    const std::complex<double> tau[],
00315    std::complex<double> C[],
00316    const int ldc,
00317    std::complex<double> work[],
00318    int* const INFO)
00319   {
00320     F77_BLAS_MANGLE(zunm2r, ZUNM2R) 
00321       (side, trans, &m, &n, &k, A, &lda, tau, C, &ldc, work, INFO);
00322   }
00323 
00324   template <>
00325   void
00326   LAPACK<int, std::complex<double> >::
00327   ORGQR (const int m,
00328    const int n,
00329    const int k,
00330    std::complex<double> A[],
00331    const int lda,
00332    std::complex<double> tau[],
00333    std::complex<double> work[],
00334    const int lwork,
00335    int* const INFO)
00336   {
00337     F77_BLAS_MANGLE(zungqr, ZUNGQR) 
00338       (&m, &n, &k, A, &lda, tau, work, &lwork, INFO);
00339   }
00340 
00341   template <>
00342   void
00343   LAPACK<int, std::complex<double> >::
00344   POTRF (const char* const uplo,
00345    const int n,
00346    std::complex<double> A[],
00347    const int lda,
00348    int* const INFO)
00349   {
00350     F77_BLAS_MANGLE(zpotrf, ZPOTRF) (uplo, &n, A, &lda, INFO);
00351   }
00352 
00353   template <>
00354   void
00355   LAPACK<int, std::complex<double> >::
00356   POTRS (const char* const uplo,
00357    const int n,
00358    const int nrhs,
00359    const std::complex<double> A[],
00360    const int lda,
00361    std::complex<double> B[],
00362    const int ldb,
00363    int* const INFO)
00364   {
00365     F77_BLAS_MANGLE(zpotrs, ZPOTRS) 
00366       (uplo, &n, &nrhs, A, &lda, B, &ldb, INFO);
00367   }
00368 
00369   template <>
00370   void
00371   LAPACK<int, std::complex<double> >::
00372   POTRI (const char* const uplo, 
00373    const int n, 
00374    std::complex<double> A[], 
00375    const int lda, 
00376    int* const INFO)
00377   {
00378     F77_BLAS_MANGLE(zpotri, ZPOTRI) (uplo, &n, A, &lda, INFO);
00379   }
00380 
00381   template <>
00382   void
00383   LAPACK<int, std::complex<double> >::
00384   LARNV (const int idist, 
00385    int iseed[],
00386    const int n,
00387    std::complex<double> x[])
00388   {
00389     F77_BLAS_MANGLE(zlarnv, ZLARNV) (&idist, iseed, &n, x);
00390   }
00391 
00392   template <>
00393   void
00394   LAPACK<int, std::complex<double> >::
00395   GESVD (const char* const jobu,
00396    const char* const jobvt,
00397    const int m,
00398    const int n,
00399    std::complex<double> A[],
00400    const int lda,
00401    double s[],
00402    std::complex<double> U[],
00403    const int ldu,
00404    std::complex<double> VT[],
00405    const int ldvt,
00406    std::complex<double> work[],
00407    const int lwork,
00408    double rwork[],
00409    int* const INFO)
00410   {
00411     F77_BLAS_MANGLE(zgesvd, ZGESVD) (jobu, jobvt, &m, &n, 
00412              A, &lda, s, 
00413              U, &ldu, VT, &ldvt, 
00414              work, &lwork, rwork, INFO);
00415   }
00416 
00417 
00418 } // namespace TSQR
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends