MechSys: linalg/lawrap.h Source File

00001 /*************************************************************************************
00002  * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
00003  * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
00004  * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
00005  *                                                                                   *
00006  * This file is part of MechSys.                                                     *
00007  *                                                                                   *
00008  * MechSys is free software; you can redistribute it and/or modify it under the      *
00009  * terms of the GNU General Public License as published by the Free Software         *
00010  * Foundation; either version 2 of the License, or (at your option) any later        *
00011  * version.                                                                          *
00012  *                                                                                   *
00013  * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
00014  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
00015  * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
00016  *                                                                                   *
00017  * You should have received a copy of the GNU General Public License along with      *
00018  * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
00019  * Fifth Floor, Boston, MA 02110-1301, USA                                           *
00020  *************************************************************************************/
00021 
00022 /****************************************************************************
00023  * Link with:                                                               *
00024  *                                                                          *
00025  *    LDFLAGS = -llapack -lf77blas   -lcblas   -latlas -lg2c                *
00026  * or                                                                       *
00027  *    LDFLAGS = -llapack -lptf77blas -lptcblas -latlas -lg2c -lpthread      *
00028  * or                                                                       *
00029  *    LDFLAGS = -llapack -lf77blas   -lcblas   -latlas -lgfortran           *
00030  * or                                                                       *
00031  *    LDFLAGS = -llapack -lptf77blas -lptcblas -latlas -lgfortran -lpthread *
00032  ****************************************************************************/
00033 
00034 #ifndef MECHSYS_LINALG_LAWRAP_H
00035 #define MECHSYS_LINALG_LAWRAP_H
00036 
00037 // STL
00038 #include <cassert>
00039 #include <map> // for Debug only
00040 
00041 // UMFPACK
00042 #ifdef USE_UMFPACK
00043   #include <umfpack.h>
00044 #endif
00045 
00046 // SCALAPACK
00047 #ifdef USE_SCALAPACK
00048   #include <mpi.h>
00049 #endif
00050 
00051 // This must be after <mpi.h>, otherwise a problem occurs.
00052 // Perhaps, mpi.h uses REAL for its own purpose
00053 // Yes, ituses:  constants.h:104:extern const Datatype REAL
00054 #ifdef HAVE_CONFIG_H
00055   #include "config.h"
00056 #else
00057   #ifndef REAL
00058     #define REAL double
00059   #endif
00060 #endif
00061 
00062 // MechSys
00063 #include "linalg/matrix.h"
00064 #include "linalg/vector.h"
00065 #include "util/exception.h"
00066 #include "util/util.h"
00067 
00068 #ifdef USE_SCALAPACK
00069   #include "linalg/scalpk.h"
00070 #endif
00071 
00072 extern "C" // {{{
00073 {
00074     // BLAS
00075     void dscal_(int const *N, double const *alpha, double *X, int const *incX);
00076 
00077     void dcopy_(int const *N, double const *X, int const *incX, double *Y, int const *incY);
00078 
00079     void daxpy_(int const *N, double const *alpha, double const *X, int const *incX, double *Y, int const *incY);
00080 
00081     double ddot_(int const *N, double const *X, int const *incX, double const *Y, int const *incY);
00082 
00083     void dsymv_(char const *Uplo, int    const *N, double const *alpha, double const *A,
00084                 int  const *lda , double const *X, int    const *incX , double const *beta,
00085                 double     *Y   , int    const *incY);
00086 
00087     void dgemv_(char   const *TransA , int    const *M   , int    const *N  , double const *alpha ,
00088                 double const *A      , int    const *lda , double const *X  , int    const *incX  ,
00089                 double const *beta   , double       *Y   , int    const *incY);
00090 
00091     void dgemm_(char   const *TransA  , char   const *TransB , int    const *M    , int    const *N   ,
00092                 int    const *K       , double const *alpha  , double const *A    , int    const *lda ,
00093                 double const *B       , int    const *ldb    , double const *beta , double       *C   , int const *ldc);
00094 
00095     void dger_(int const *M   , int    const *N, double const *alpha, double const *X,
00096                int const *incX, double const *Y, int    const *incY , double       *A, int const *lda);
00097 
00098     // DGESV
00099 #if defined (USE_LAPACK) || defined (USE_GOTOBLAS)
00100   #ifdef USE_FLOAT
00101     void sgesv_(int *Np, int *NRHSp,
00102                 float *A, int *LDAp, int *IPIVp,
00103                 float *B, int *LDBp, int *INFOp);
00104   #else
00105     void dgesv_(int *Np, int *NRHSp,
00106                 double *A, int *LDAp, int *IPIVp,
00107                 double *B, int *LDBp, int *INFOp);
00108   #endif // #ifdef USE_FLOAT
00109 #endif // #if defined (USE_LAPACK) || defined (USE_GOTOBLAS)
00110 
00111     // DGTSV and EIGENPROBLEMS
00112 #ifdef USE_LAPACK
00113   #ifdef USE_FLOAT
00114     float slamch_ (char *CMACHp);
00115     void  sgtsv_  (int *Np, int *NRHSp, float *DL, float *D, float *DU, float *B, int *LDBp, int *INFOp);
00116     void  ssyevr_ (char *JOBZp, char *RANGEp, char *UPLOp, int *Np, float *A, int *LDAp, float *VLp, float *VUp,
00117                    int *ILp, int *IUp, float *ABSTOLp, int *Mp, float *W, float *Z, int *LDZp, int *ISUPPZ,
00118                    float *WORK, int *LWORKp, int *IWORK, int *LIWORKp, int *INFOp);
00119   #else
00120     double dlamch_ (char *CMACHp);
00121     void   dgtsv_  (int *Np, int *NRHSp, double *DL, double *D, double *DU, double *B, int *LDBp, int *INFOp);
00122     void   dsyevr_ (char *JOBZp, char *RANGEp, char *UPLOp, int *Np, double *A, int *LDAp, double *VLp, double *VUp,
00123                     int *ILp, int *IUp, double *ABSTOLp, int *Mp, double *W, double *Z, int *LDZp, int *ISUPPZ,
00124                     double *WORK, int *LWORKp, int *IWORK, int *LIWORKp, int *INFOp);
00125   #endif // #ifdef USE_FLOAT
00126 #endif // #ifdef USE_LAPACK
00127 
00128 } // }}}
00129 
00130 namespace LinAlg
00131 {
00132 
00194 // {{{ BLAS
00200 // ######################################################################################################################## BLAS
00201 
00202 // {{{ Level1 (vector-vector)
00208 // =============================================================================== vector-vector (Level 1)
00209 
00210 // -------------------------------------------------------------------------- Scal
00212 inline void Scal(REAL   const   a, 
00213                  Vector<REAL> & X) 
00214 {
00215     int N = X.Size();
00216 #ifdef USE_FLOAT
00217     throw new Fatal(_("LinAlg::Scal not yet."));
00218     //cblas_sscal(N,a,X.GetPtr(),1);
00219 #else
00220     int incX = 1;
00221     dscal_(&N,&a,X.GetPtr(),&incX);
00222 #endif // #ifdef USE_FLOAT
00223 }
00224 
00225 // -------------------------------------------------------------------------- Copy
00227 inline void Copy(Vector<REAL> const & X, 
00228                  Vector<REAL>       & Y) 
00229 {
00230     int N = X.Size();
00231     assert(Y.Size()==N);
00232 #ifdef USE_FLOAT
00233     throw new Fatal(_("LinAlg::Copy not yet."));
00234     //cblas_scopy(N,X.GetPtr(),1,Y.GetPtr(),1);
00235 #else
00236     int incX = 1;
00237     int incY = 1;
00238     dcopy_(&N,X.GetPtr(),&incX,Y.GetPtr(),&incY);
00239 #endif // #ifdef USE_FLOAT
00240 }
00241 
00242 // -------------------------------------------------------------------------- Axpy
00244 inline void Axpy(REAL         const   a, 
00245                  Vector<REAL> const & X, 
00246                  Vector<REAL>       & Y) 
00247 {
00248     int N = X.Size();
00249     assert(Y.Size()==N);
00250 #ifdef USE_FLOAT
00251     throw new Fatal(_("LinAlg::Axpy not yet."));
00252     //cblas_saxpy(N,a,X.GetPtr(),1,Y.GetPtr(),1);
00253 #else
00254     int incX = 1;
00255     int incY = 1;
00256     daxpy_(&N,&a,X.GetPtr(),&incX,Y.GetPtr(),&incY);
00257 #endif // #ifdef USE_FLOAT
00258 }
00259 
00260 // --------------------------------------------------------------------------- Dot
00262 
00263 #ifdef USE_FLOAT
00264 inline float  Dot(Vector<REAL> const & X, 
00265                   Vector<REAL> const & Y) 
00266 #else
00267 inline double Dot(Vector<REAL> const & X, 
00268                   Vector<REAL> const & Y) 
00269 #endif // #ifdef USE_FLOAT
00270 {
00271     int N = X.Size();
00272     assert(Y.Size()==N);
00273 #ifdef USE_FLOAT
00274     throw new Fatal(_("LinAlg::Dot not yet."));
00275     //return cblas_sdot(N,X.GetPtr(),1,Y.GetPtr(),1);
00276 #else
00277     int incX = 1;
00278     int incY = 1;
00279     return ddot_(&N,X.GetPtr(),&incX,Y.GetPtr(),&incY);
00280 #endif // #ifdef USE_FLOAT
00281 }
00282 
00283 // -------------------------------------------------------------------------- Norm
00285 inline REAL Norm(Vector<REAL> const & X) 
00286 {
00287     return sqrt(LinAlg::Dot(X,X));
00288 }
00289 
00290 // ---------------------------------------------------------------------- CopyScal
00292 inline void CopyScal(REAL         const   a, 
00293                      Vector<REAL> const & X, 
00294                      Vector<REAL>       & Y) 
00295 {
00296     assert(Y.Size()==X.Size());
00297     Copy(X,Y); // Y <- X;
00298     Scal(a,Y); // Y <- aY = aX
00299 }
00300 
00301 // ---------------------------------------------------------------------- AddScaled
00303 inline void AddScaled(REAL         const   a,
00304                       Vector<REAL> const & X,
00305                       REAL         const   b,
00306                       Vector<REAL> const & Y,
00307                       Vector<REAL>       & Z)
00308 {
00309     assert(Y.Size()==X.Size());
00310     Z.SetValues(0.0); // Z <- 0.0
00311     Axpy(a,X,Z);      // Z <- a*X + 0.0
00312     Axpy(b,Y,Z);      // Z <- b*Y + a*X
00313 }
00314  // end of sub-group Level1
00316 // }}}
00317 
00318 // {{{ Level2 (matrix-vector)
00324 // =============================================================================== matrix-vector (Level 2)
00325 
00326 // -------------------------------------------------------------------------- Symv
00328 
00343 inline void Symv(REAL          const   a, 
00344                  Matrix<REAL>  const & A, 
00345                  Vector<REAL>  const & X, 
00346                  REAL          const   b, 
00347                  Vector<REAL>        & Y) 
00348 {
00349     int N = A.Rows();
00350     assert(A.Cols()==N);
00351     assert(X.Size()==N);
00352     assert(Y.Size()==N);
00353 #ifdef USE_FLOAT
00354     throw new Fatal(_("LinAlg::Symv not yet."));
00355     //cblas_ssymv(CblasColMajor, // CBLAS_ORDER
00356                 //CblasUpper,    // CBLAS_UPLO
00357                 //N,             // N
00358                 //a,             // Alpha
00359                 //A.GetPtr(),    // const float  * A
00360                 //N,             // LDA
00361                 //X.GetPtr(),    // const float  * X
00362                 //1,             // incX
00363                 //b,             // Beta
00364                 //Y.GetPtr(),    // float  * Y
00365                 //1);            // incY
00366 #else
00367     char Uplo = 'U';
00368     int  incX = 1;
00369     int  incY = 1;
00370     dsymv_(&Uplo,         // CBLAS_UPLO
00371            &N,            // N
00372            &a,            // Alpha
00373            A.GetPtr(),    // const double * A
00374            &N,            // LDA
00375            X.GetPtr(),    // const double * X
00376            &incX,         // incX
00377            &b,            // Beta
00378            Y.GetPtr(),    // double * Y
00379            &incY);        // incY
00380 #endif // #ifdef USE_FLOAT
00381 }
00382 
00383 // -------------------------------------------------------------------------- Gemv
00385 
00401 inline void Gemv(REAL          const   a, 
00402                  Matrix<REAL>  const & A, 
00403                  Vector<REAL>  const & X, 
00404                  REAL          const   b, 
00405                  Vector<REAL>        & Y) 
00406 {
00407     int M = A.Rows();
00408     int N = A.Cols();
00409     assert(X.Size()==N);
00410     assert(Y.Size()==M);
00411 #ifdef USE_FLOAT
00412     throw new Fatal(_("LinAlg::Gemv not yet."));
00413     //cblas_sgemv(CblasColMajor, // CBLAS_ORDER
00414                 //CblasNoTrans,  // CBLAS_TRANSPOSE A
00415                 //M,             // M
00416                 //N,             // N
00417                 //a,             // Alpha
00418                 //A.GetPtr(),    // const float  * A
00419                 //M,             // LDA
00420                 //X.GetPtr(),    // const float  * X
00421                 //1,             // incX
00422                 //b,             // Beta
00423                 //Y.GetPtr(),    // float  * Y
00424                 //1);            // incY
00425 #else
00426     char TransA = 'N';
00427     int  incX   = 1;
00428     int  incY   = 1;
00429     dgemv_(&TransA,        // CBLAS_TRANSPOSE A
00430            &M,             // M
00431            &N,             // N
00432            &a,             // Alpha
00433            A.GetPtr(),     // const double * A
00434            &M,             // LDA
00435            X.GetPtr(),     // const double * X
00436            &incX,          // incX
00437            &b,             // Beta
00438            Y.GetPtr(),     // double * Y
00439            &incY);         // incY
00440 #endif // #ifdef USE_FLOAT
00441 }
00442 
00443 // ------------------------------------------------------------------------- Gemtv
00445 inline void Gemtv(REAL          const   a, 
00446                   Matrix<REAL>  const & A, 
00447                   Vector<REAL>  const & X, 
00448                   REAL          const   b, 
00449                   Vector<REAL>        & Y) 
00450 {
00451     int M = A.Rows();
00452     int N = A.Cols();
00453     assert(X.Size()==N);
00454     assert(Y.Size()==M);
00455 #ifdef USE_FLOAT
00456     throw new Fatal(_("LinAlg::Gemtv not yet."));
00457     //cblas_sgemv(CblasColMajor, // CBLAS_ORDER
00458                 //CblasTrans,    // CBLAS_TRANSPOSE A
00459                 //M,             // M
00460                 //N,             // N
00461                 //a,             // Alpha
00462                 //A.GetPtr(),    // const float  * A
00463                 //M,             // LDA
00464                 //X.GetPtr(),    // const float  * X
00465                 //1,             // incX
00466                 //b,             // Beta
00467                 //Y.GetPtr(),    // float  * Y
00468                 //1);            // incY
00469 #else
00470     char TransA = 'T';
00471     int  incX   = 1;
00472     int  incY   = 1;
00473     dgemv_(&TransA,        // CBLAS_TRANSPOSE A
00474            &M,             // M
00475            &N,             // N
00476            &a,             // Alpha
00477            A.GetPtr(),     // const double * A
00478            &M,             // LDA
00479            X.GetPtr(),     // const double * X
00480            &incX,          // incX
00481            &b,             // Beta
00482            Y.GetPtr(),     // double * Y
00483            &incY);         // incY
00484 #endif // #ifdef USE_FLOAT
00485 }
00486 
00487 // --------------------------------------------------------------------------- Ger
00489 
00502 inline void Ger(REAL         const   a, 
00503                 Vector<REAL> const & X, 
00504                 Vector<REAL> const & Y, 
00505                 Matrix<REAL>       & A) 
00506 {
00507     int M = A.Rows();
00508     int N = A.Cols();
00509     assert(X.Size()==M);
00510     assert(Y.Size()==N);
00511 #ifdef USE_FLOAT
00512     throw new Fatal(_("LinAlg::Ger not yet."));
00513     //cblas_sger(CblasColMajor, // CBLAS_ORDER
00514                //M,             // M
00515                //N,             // N
00516                //a,             // Alpha
00517                //X.GetPtr(),    // const float  * X
00518                //1,             // incX
00519                //Y.GetPtr(),    // const float  * Y
00520                //1,             // incY
00521                //A.GetPtr(),    // float  * A
00522                //M);            // LDA
00523 #else
00524     int incX = 1;
00525     int incY = 1;
00526     dger_(&M,             // M
00527           &N,             // N
00528           &a,             // Alpha
00529           X.GetPtr(),     // const double * X
00530           &incX,          // incX
00531           Y.GetPtr(),     // const double * Y
00532           &incY,          // incY
00533           A.GetPtr(),     // double * A
00534           &M);            // LDA
00535 #endif // #ifdef USE_FLOAT
00536 }
00537 
00538 // ========================================================================== Derived
00539 
00540 // -------------------------------------------------------------------------- Gemvpmv
00542 
00549 inline void Gemvpmv(REAL          const   a, 
00550                     Matrix<REAL>  const & A, 
00551                     Vector<REAL>  const & X, 
00552                     REAL          const   b, 
00553                     Matrix<REAL>  const & B, 
00554                     Vector<REAL>  const & Y, 
00555                     Vector<REAL>        & Z) 
00556 {
00557 #ifndef NDEBUG
00558     int M = A.Rows();
00559     int N = A.Cols();
00560     int P = B.Cols();
00561     assert(B.Rows()==M);
00562     assert(X.Size()==N);
00563     assert(Y.Size()==P);
00564     assert(Z.Size()==M);
00565 #endif
00566     Gemv(a,A,X, 0.0,Z);  // Z <- a*A*X + 0*Z
00567     Gemv(b,B,Y, 1.0,Z);  // Z <- b*B*Y + 1*Z
00568 }
00569 
00570 // ============ LAPACK
00571  // end of sub-group Level2
00573 // }}}
00574 
00575 // {{{ Level3 (matrix-matrix)
00581 // =============================================================================== matrix-matrix (Level 3)
00582 
00583 // -------------------------------------------------------------------------- Gemm
00585 
00599 inline void Gemm(REAL         const   a, 
00600                  Matrix<REAL> const & A, 
00601                  Matrix<REAL> const & B, 
00602                  REAL         const   b, 
00603                  Matrix<REAL>       & C) 
00604 {
00605     int M = A.Rows();
00606     int K = A.Cols();
00607     int N = B.Cols();
00608     assert(B.Rows()==K);
00609     assert(C.Rows()==M); assert(C.Cols()==N);
00610 #ifdef USE_FLOAT
00611     throw new Fatal(_("LinAlg::Gemm not yet."));
00612     //cblas_sgemm(CblasColMajor, // CBLAS_ORDER
00613                 //CblasNoTrans,  // CBLAS_TRANSPOSE A
00614                 //CblasNoTrans,  // CBLAS_TRANSPOSE B
00615                 //M,             // M
00616                 //N,             // N
00617                 //K,             // K
00618                 //a,             // Alpha
00619                 //A.GetPtr(),    // const float  * A
00620                 //M,             // LDA
00621                 //B.GetPtr(),    // const float  * B
00622                 //K,             // LDB
00623                 //b,             // Beta
00624                 //C.GetPtr(),    // float  * C
00625                 //M);            // LDC
00626 #else
00627     char TransA='N';
00628     char TransB='N';
00629     dgemm_(&TransA,       // CBLAS_TRANSPOSE A
00630            &TransB,       // CBLAS_TRANSPOSE B
00631            &M,            // M
00632            &N,            // N
00633            &K,            // K
00634            &a,            // Alpha
00635            A.GetPtr(),    // const double * A
00636            &M,            // LDA
00637            B.GetPtr(),    // const double * B
00638            &K,            // LDB
00639            &b,            // Beta
00640            C.GetPtr(),    // double * C
00641            &M);           // LDC
00642 #endif // #ifdef USE_FLOAT
00643 }
00644 
00645 // ------------------------------------------------------------------------- Gemtm
00647 inline void Gemtm(REAL         const   a, 
00648                   Matrix<REAL> const & A, 
00649                   Matrix<REAL> const & B, 
00650                   REAL         const   b, 
00651                   Matrix<REAL>       & C) 
00652 {
00653     int K = A.Rows();
00654     int M = A.Cols();
00655     int N = B.Cols();
00656     assert(B.Rows()==K);
00657     assert(C.Rows()==M); assert(C.Cols()==N);
00658 #ifdef USE_FLOAT
00659     throw new Fatal(_("LinAlg::Gemtm not yet."));
00660     //cblas_sgemm(CblasColMajor, // CBLAS_ORDER
00661                 //CblasTrans,    // CBLAS_TRANSPOSE A
00662                 //CblasNoTrans,  // CBLAS_TRANSPOSE B
00663                 //M,             // M
00664                 //N,             // N
00665                 //K,             // K
00666                 //a,             // Alpha
00667                 //A.GetPtr(),    // const float  * A
00668                 //K,             // LDA
00669                 //B.GetPtr(),    // const float  * B
00670                 //K,             // LDB
00671                 //b,             // Beta
00672                 //C.GetPtr(),    // float  * C
00673                 //M);            // LDC
00674 #else
00675     char TransA='T';
00676     char TransB='N';
00677     dgemm_(&TransA,       // CBLAS_TRANSPOSE A
00678            &TransB,       // CBLAS_TRANSPOSE B
00679            &M,            // M
00680            &N,            // N
00681            &K,            // K
00682            &a,            // Alpha
00683            A.GetPtr(),    // const double * A
00684            &K,            // LDA
00685            B.GetPtr(),    // const double * B
00686            &K,            // LDB
00687            &b,            // Beta
00688            C.GetPtr(),    // double * C
00689            &M);           // LDC
00690 #endif // #ifdef USE_FLOAT
00691 }
00692 
00693 // ------------------------------------------------------------------------- Gemmt
00695 inline void Gemmt(REAL         const   a, 
00696                   Matrix<REAL> const & A, 
00697                   Matrix<REAL> const & B, 
00698                   REAL         const   b, 
00699                   Matrix<REAL>       & C) 
00700 {
00701     int M = A.Rows();
00702     int K = A.Cols();
00703     int N = B.Rows();
00704     assert(B.Cols()==K);
00705     assert(C.Rows()==M); assert(C.Cols()==N);
00706 #ifdef USE_FLOAT
00707     throw new Fatal(_("LinAlg::Gemmt not yet."));
00708     //cblas_sgemm(CblasColMajor, // CBLAS_ORDER
00709                 //CblasNoTrans,  // CBLAS_TRANSPOSE A
00710                 //CblasTrans,    // CBLAS_TRANSPOSE B
00711                 //M,             // M
00712                 //N,             // N
00713                 //K,             // K
00714                 //a,             // Alpha
00715                 //A.GetPtr(),    // const float  * A
00716                 //M,             // LDA
00717                 //B.GetPtr(),    // const float  * B
00718                 //N,             // LDB
00719                 //b,             // Beta
00720                 //C.GetPtr(),    // float  * C
00721                 //M);            // LDC
00722 #else
00723     char TransA='N';
00724     char TransB='T';
00725     dgemm_(&TransA,       // CBLAS_TRANSPOSE A
00726            &TransB,       // CBLAS_TRANSPOSE B
00727            &M,            // M
00728            &N,            // N
00729            &K,            // K
00730            &a,            // Alpha
00731            A.GetPtr(),    // const double * A
00732            &M,            // LDA
00733            B.GetPtr(),    // const double * B
00734            &N,            // LDB
00735            &b,            // Beta
00736            C.GetPtr(),    // double * C
00737            &M);           // LDC
00738 #endif // #ifdef USE_FLOAT
00739 }
00740 
00741 // ------------------------------------------------------------------------ Gemtmt
00743 inline void Gemtmt(REAL         const   a, 
00744                    Matrix<REAL> const & A, 
00745                    Matrix<REAL> const & B, 
00746                    REAL         const   b, 
00747                    Matrix<REAL>       & C) 
00748 {
00749     int K = A.Rows();
00750     int M = A.Cols();
00751     int N = B.Rows();
00752     assert(B.Cols()==K);
00753     assert(C.Rows()==M); assert(C.Cols()==N);
00754 #ifdef USE_FLOAT
00755     throw new Fatal(_("LinAlg::Gemtmt not yet."));
00756     //cblas_sgemm(CblasColMajor, // CBLAS_ORDER
00757                 //CblasTrans,    // CBLAS_TRANSPOSE A
00758                 //CblasTrans,    // CBLAS_TRANSPOSE B
00759                 //M,             // M
00760                 //N,             // N
00761                 //K,             // K
00762                 //a,             // Alpha
00763                 //A.GetPtr(),    // const float  * A
00764                 //K,             // LDA
00765                 //B.GetPtr(),    // const float  * B
00766                 //N,             // LDB
00767                 //b,             // Beta
00768                 //C.GetPtr(),    // float  * C
00769                 //M);            // LDC
00770 #else
00771     char TransA='T';
00772     char TransB='T';
00773     dgemm_(&TransA,       // CBLAS_TRANSPOSE A
00774            &TransB,       // CBLAS_TRANSPOSE B
00775            &M,            // M
00776            &N,            // N
00777            &K,            // K
00778            &a,            // Alpha
00779            A.GetPtr(),    // const double * A
00780            &K,            // LDA
00781            B.GetPtr(),    // const double * B
00782            &N,            // LDB
00783            &b,            // Beta
00784            C.GetPtr(),    // double * C
00785            &M);           // LDC
00786 #endif // #ifdef USE_FLOAT
00787 }
00788  // end of sub-group Level3
00790 // }}}
00791  // end of group BLAS }}}
00793 
00794 
00795 // ========================================================================= Solver
00796 
00797 // -------------------------------------------------------------------------- Gesv
00799 
00808 inline int Gesv(Matrix<REAL> & A,                 
00809                 Vector<REAL> & Y,                 
00810                 bool           DoPreserveA=false) 
00811 {
00812     int info = 0;
00813 
00814 #if (defined (USE_LAPACK) || defined (USE_GOTOBLAS))
00815 
00816     // {{{
00817     
00818     // Declare a pointer to a temporary A matrix
00819     Matrix<REAL> * tmp_A; // tmp_A <- A
00820 
00821     // Assign/allocate tmp_A
00822     if (DoPreserveA==true) tmp_A = new Matrix<REAL> (A); // tmp_A <- A
00823     else                   tmp_A = &A;
00824 
00825     // Solve
00826     int    N = tmp_A->Rows();
00827     int NRHS = 1;        //int NRHS = Y.Cols();
00828     assert(tmp_A->Cols()==N);
00829     assert(Y.Size()==N); //assert(Y.Rows()==N);
00830     int * ipiv = new int [N];
00831 
00832   #ifdef USE_FLOAT
00833     sgesv_(&N,                 // A(N,N)
00834            &NRHS,              // {Y}(N,NRHS) (RHS: Right Hand Side) 
00835            tmp_A->GetPtr(),    // double * A
00836            &N,                 // LDA
00837            ipiv,               // Pivot Indices
00838            Y.GetPtr(),         // double * Y
00839            &N,                 // LDY
00840            &info);
00841   #else
00842     dgesv_(&N,                 // A(N,N)
00843            &NRHS,              // {Y}(N,NRHS) (RHS: Right Hand Side) 
00844            tmp_A->GetPtr(),    // double * A
00845            &N,                 // LDA
00846            ipiv,               // Pivot Indices
00847            Y.GetPtr(),         // double * Y
00848            &N,                 // LDY
00849            &info);
00850   #endif // #ifdef USE_FLOAT
00851 
00852     // Clean up
00853     delete [] ipiv;
00854     if (DoPreserveA==true) delete tmp_A;
00855 
00856     // }}}
00857     
00858 #elif (defined (USE_UMFPACK))
00859 
00860     // {{{
00861 
00862     // This part will preserve A matrix, because of temporary copies
00863 
00864     // Tolerance
00865     const double ZEROELEM = 1.0e-9;
00866 
00867     // Size
00868     int m = A.Rows();
00869     int n = A.Cols();
00870 
00871     // Number of nonzero elements
00872     int num_nonzero=0;
00873     for (int j=0; j<n; ++j)
00874     for (int i=0; i<m; ++i)
00875         if (fabs(A(i,j))>ZEROELEM) num_nonzero++;
00876 
00877     // Allocate arrays
00878     int    * Ap = new int    [n+1];
00879     int    * Ai = new int    [num_nonzero];
00880     double * Ax = new double [num_nonzero];
00881     for (int j=0; j<n; ++j) Ap[j]=-1; // <<< not set yet
00882 
00883     // Fill arrays
00884     int k=0;
00885     for (int j=0; j<n; ++j)
00886     for (int i=0; i<m; ++i)
00887     {
00888         if (fabs(A(i,j))>ZEROELEM)
00889         {
00890             if (Ap[j]==-1) // not set yet
00891                 Ap[j]=k;
00892             Ai[k] = i;
00893             Ax[k] = A(i,j);
00894             k++;
00895         }
00896     }
00897     Ap[n] = num_nonzero;
00898 
00899     #ifdef UMFPACK_DEBUG
00900     // Debug
00901     cout << "num_nonzero = " << num_nonzero << endl;
00902     cout << "Ap = "; for (int j=0; j<n+1; ++j) cout<<Ap[j]<<"  "; cout << endl;
00903     cout << "Ai = "; for (int i=0; i<k;   ++i) cout<<Ai[i]<<"  "; cout << endl;
00904     cout << "Ax = "; for (int i=0; i<k;   ++i) cout<<Ax[i]<<"  "; cout << endl;
00905     #endif
00906 
00907     // Solve
00908     Vector<REAL> b(Y);
00909     double *null = (double *)NULL;
00910     void *Symbolic, *Numeric;
00911     info  = umfpack_di_symbolic      (n, n, Ap, Ai, Ax, &Symbolic, null, null);
00912     info += umfpack_di_numeric       (Ap, Ai, Ax, Symbolic, &Numeric, null, null);
00913             umfpack_di_free_symbolic (&Symbolic);
00914     info += umfpack_di_solve         (UMFPACK_A, Ap, Ai, Ax, Y.GetPtr(), b.GetPtr(), Numeric, null, null);
00915             umfpack_di_free_numeric  (&Numeric);
00916 
00917     // Clean up
00918     delete [] Ap;
00919     delete [] Ai;
00920     delete [] Ax;
00921 
00922     // }}}
00923 
00924 #elif (defined (USE_SCALAPACK))
00925 
00926     // {{{
00927     
00928     // This part will preserve A matrix, because of temporary copies
00929     
00930     // Information
00931     int size = MPI::COMM_WORLD.Get_size(); // Number of process
00932 
00933     // Factor the number of processors into a rectangle quadratic shape preferred
00934     int nprow = -1;
00935     int npcol = -1;
00936     Util::FindBestSquare(size, nprow,npcol);
00937     if ((nprow<1) || (npcol<1))
00938     {
00939         MPI::Finalize();
00940         throw new Fatal(_("LinAlg::Gesv::ScaLAPACK: Could not find nprow=%d and npcol=%d"),nprow,npcol);
00941     }
00942 
00943     // Determine a good block size
00944     int N    = A.Rows();
00945     int NRHS = 1;
00946     int nb   = 2;
00947     if ((N/nprow)<(N/npcol)) nb=(N/nprow);
00948     else                     nb=(N/npcol);
00949     if (nb<1)  nb=1;
00950     if (nb>64) nb=64;
00951 
00952     // Debug
00953     #ifdef SCALAPACK_DEBUG
00954     nprow=2;
00955     npcol=3;
00956     nb   =2;
00957     if (size!=6)
00958     {
00959         MPI::Finalize();
00960         throw new Fatal(_("LinAlg::Gesv::ScaLAPACK: The number of process for DEBUG must be equal to 6"));
00961     }
00962     #endif
00963 
00964     // Check
00965     if (nb>N) nb=N;
00966     if (nprow*npcol>size)
00967     {
00968         MPI::Finalize();
00969         throw new Fatal(_("LinAlg::Gesv::ScaLAPACK: There is not enough processes (%d<%d) to make a %d-by-%d process grid"),size,nprow*npcol,nprow,npcol);
00970     }
00971 
00972     // Initialize BLACS process grid
00973     int iam,   nprocs, ictxt;
00974     int myrow, mycol;
00975     Cblacs_pinfo    (&iam, &nprocs);
00976     Cblacs_get      (-1, 0, &ictxt);
00977     Cblacs_gridinit (&ictxt, /*order*/"Col-major", nprow, npcol);
00978     Cblacs_gridinfo (ictxt, &nprow, &npcol, &myrow, &mycol);
00979 
00980     #ifdef SCALAPACK_DEBUG
00981     // Debug
00982     String msg, buf;
00983     msg.Printf("Hi, I'am {%d,%d}: nprow=%d, npcol=%d, nb=%d",myrow,mycol,nprow,npcol,nb);
00984     std::map<int,String> letters;
00985     letters[19]="S";  letters[-19]="-S";
00986     letters[ 3]="C";  letters[ -3]="-C";
00987     letters[ 1]="A";  letters[ -1]="-A";
00988     letters[12]="L";  letters[-12]="-L";
00989     letters[16]="P";  letters[-16]="-P";
00990     letters[11]="K";  letters[-11]="-K";
00991     msg.append(" A =\n");
00992     int rank = MPI::COMM_WORLD.Get_rank();
00993     if (rank==0)
00994     {
00995         for (int i=0; i<A.Rows(); ++i)
00996         {
00997             for (int j=0; j<A.Rows(); ++j)
00998             {
00999                 buf.Printf("   %2s",letters[static_cast<int>(A(i,j))].GetSTL().c_str());
01000                 msg.append(buf);
01001             }
01002             msg.append("\n");
01003         }
01004     }
01005     #endif
01006     
01007     // Work only the process in the process grid
01008     if ((myrow>-1)&&(mycol>-1)&&(myrow<nprow)&&(mycol<npcol))
01009     {
01010         // Compute the size of the local matrices
01011         int isrc    = 0;
01012         int m_loc_A = numroc_(&N   , &nb, &myrow, &isrc, &nprow); // No of rows in local data block
01013         int n_loc_A = numroc_(&N   , &nb, &mycol, &isrc, &npcol); // No of columns in local data block
01014         int n_loc_b = numroc_(&NRHS, &nb, &mycol, &isrc, &npcol); // No of columns in local data block
01015 
01016         // m_loc_b must be equal to m_loc_A (TODO: check this ?)
01017         
01018         #ifdef SCALAPACK_DEBUG
01019         // Debug
01020         buf.Printf(" m_loc_A=%d, n_loc_A=%d, n_loc_b=%d\n", m_loc_A,n_loc_A,n_loc_b);
01021         msg.append(buf);
01022         #endif
01023 
01024         // Allocate and fill the matrices [loc_A] and [loc_b]
01025         double * loc_A = new double [m_loc_A*n_loc_A]; // local A (sub) array
01026         double * loc_b = NULL;
01027         if (n_loc_b>0)
01028             loc_b = new double [m_loc_A*n_loc_b]; // local b (sub) array
01029 
01030         // Initialize the array descriptor for A and b arrays
01031         int loc_ld_A = Util::Max(1,m_loc_A); // local leading dimension
01032         int loc_ld_b = loc_ld_A;
01033         int desc_A[9];
01034         int desc_b[9];
01035         descinit_(desc_A, &N, &N   , &nb, &nb, &isrc, &isrc, &ictxt, &loc_ld_A, &info);
01036         descinit_(desc_b, &N, &NRHS, &nb, &nb, &isrc, &isrc, &ictxt, &loc_ld_b, &info);
01037 
01038         // Fill [loc_A] and [loc_b] (sub) arrays
01039         for (int i_loc=0; i_loc<m_loc_A; i_loc++)
01040         {
01041             int i_glob = IndxL2G(i_loc, nb, myrow, isrc, nprow);
01042             for (int j_loc=0; j_loc<n_loc_A; j_loc++)
01043             {
01044                 int j_glob = IndxL2G(j_loc, nb, mycol, isrc, npcol);
01045                 #ifdef SCALAPACK_DEBUG
01046                 // Debug
01047                 buf.Printf("  [%d,%d]=A(%d,%d)=%2s",i_loc,j_loc,i_glob,j_glob,letters[static_cast<int>(A(i_glob,j_glob))].GetSTL().c_str());
01048                 msg.append(buf);
01049                 #endif
01050                 // Fill [loc_A]
01051                 loc_A[j_loc*m_loc_A+i_loc] = A(i_glob,j_glob);
01052             }
01053             #ifdef SCALAPACK_DEBUG
01054             // Debug
01055             msg.append("\n");
01056             #endif
01057             // Fill [loc_b]
01058             if (n_loc_b>0) loc_b[i_loc] = Y(i_glob);
01059         }
01060 
01061         // Allocate a 'pivot' array
01062         int * ippiv = new int [m_loc_A+nb];
01063 
01064         // Call ScaLAPACK PDGESV routine
01065         int ione = 1;
01066         #ifdef SCALAPACK_DEBUG
01067         double tini = MPI::Wtime();
01068         #endif
01069 
01070         pdgesv_(&N, &NRHS, loc_A, &ione, &ione, desc_A, ippiv, loc_b, &ione, &ione, desc_b, &info);
01071 
01072         #ifdef SCALAPACK_DEBUG
01073         // Debug
01074         double tfin = MPI::Wtime();
01075         if (n_loc_b>0)
01076         {
01077             msg.append("\n loc_b = ");
01078             for (int i_loc=0; i_loc<m_loc_A; i_loc++)
01079             {
01080                 buf.Printf("%12f ",loc_b[i_loc]);
01081                 msg.append(buf);
01082             }
01083             msg.append("\n\n");
01084         }
01085         buf.Printf(" Elapsed time = %f\n",tfin-tini);
01086         msg.append(buf);
01087         #endif
01088 
01089         // Fill global solution vector (all-to-gather)
01090         int lld_Y = N;
01091         int desc_Y[9]; 
01092         int my = N*nprow;
01093         int ny = NRHS*npcol;
01094         descinit_(desc_Y, &my, &ny, &N, &NRHS, &isrc, &isrc, &ictxt, &lld_Y, &info);
01095         for (int i=0; i<nprow; ++i)
01096         {
01097             int istart = i*N+1;
01098             for (int j=0; j<npcol; ++j)
01099             {
01100                 int jstart = j*NRHS+1;
01101                 Cpdgemr2d(N,NRHS,loc_b,1,1,desc_b, Y.GetPtr(),istart,jstart,desc_Y,ictxt);
01102             }
01103         }
01104 
01105         // Wait until Y is filled
01106         Cblacs_barrier(ictxt, "All");
01107 
01108         #ifdef SCALAPACK_DEBUG
01109         // Debug
01110         msg.append("\n\n Solution: Y = ");
01111         for (int i=0; i<Y.Size(); ++i)
01112         {
01113             buf.Printf("%e ",Y(i));
01114             msg.append(buf);
01115         }
01116         msg.append("\n\n");
01117         std::cout<<msg<<std::endl;
01118         #endif
01119 
01120         // Clean up
01121         delete [] loc_A;
01122         if (n_loc_b>0)
01123             delete [] loc_b;
01124         delete [] ippiv;
01125 
01126         // Grid exit
01127         Cblacs_gridexit(0);
01128     }
01129 
01130     // }}}
01131 
01132 #else
01133     throw new Fatal(_("LinAlg::Gesv: Low level solver {dgesv_, pdgesv_, ...} is not available"));
01134 #endif
01135 
01136     if (info!=0)
01137     {
01138 #ifdef USE_SCALAPACK
01139         MPI::Finalize();
01140 #endif
01141         throw new Fatal(_("LinAlg::Gesv: info=%d means that dgesv routine failed"),info);
01142     }
01143     return info;
01144 }
01145 
01146 // ============================================================= Tridiagonal solver
01147 
01148 #ifdef USE_LAPACK
01149 
01150 // {{{
01151 
01152 // -------------------------------------------------------------------------- Gtsv
01154 
01173 inline int Gtsv(Vector<REAL> & DL, 
01174                 Vector<REAL> & D,  
01175                 Vector<REAL> & DU, 
01176                 Vector<REAL> & Y)  
01177 {
01178     int    N = D.Size();
01179     int NRHS = 1;          //int NRHS = Y.Cols();
01180     assert(Y.Size()==N);   //assert(Y.Rows()==N);
01181     assert(DL.Size()==N-1);
01182     assert(DU.Size()==N-1);
01183     int info;
01184 #ifdef USE_FLOAT
01185     sgtsv_(&N,          // A(N,N)
01186            &NRHS,       // {Y}(N,NRHS) (RHS: Right Hand Side)
01187            DL.GetPtr(), // float  * DL
01188            D.GetPtr(),  // float  * D
01189            DU.GetPtr(), // float  * DU
01190            Y.GetPtr(),  // float  * Y
01191            &N,          // LDY
01192            &info);      // int * INFO
01193 #else
01194     dgtsv_(&N,          // A(N,N)
01195            &NRHS,       // {Y}(N,NRHS) (RHS: Right Hand Side)
01196            DL.GetPtr(), // double * DL
01197            D.GetPtr(),  // double * D
01198            DU.GetPtr(), // double * DU
01199            Y.GetPtr(),  // double * Y
01200            &N,          // LDY
01201            &info);      // int * INFO
01202 #endif // #ifdef USE_FLOAT
01203     return info;
01204 }
01205 
01206 // }}}
01207 
01208 #endif
01209 
01210 // ================================================================== Eigenproblems
01211 
01212 #ifdef USE_LAPACK
01213 
01214 // {{{
01215 
01216 // ------------------------------------------------------------------------- Syevr
01218 inline int Syevr(Matrix<REAL> & A, 
01219                  Vector<REAL> & L, 
01220                  Matrix<REAL> & Q) 
01221 {
01222     int N = A.Rows();
01223     assert(A.Cols()==N);
01224     assert(L.Size()==N);
01225     assert(Q.Rows()==N); assert(Q.Cols()==N);
01226     char     jobz = 'V';
01227     char     range = 'A';
01228     char     uplo = 'L';
01229     REAL     vl,vu;
01230     int      il,iu;
01231     char     cmach = 'S';
01232 #ifdef USE_FLOAT
01233     REAL     abstol = slamch_(&cmach);
01234 #else
01235     REAL     abstol = dlamch_(&cmach);
01236 #endif
01237     int      m;
01238     int    * isuppz = new int [2*N];
01239     int      lwork = 26*N;
01240     REAL   * work = new REAL [lwork];
01241     int      liwork = 10*N;
01242     int    * iwork = new int [liwork];
01243     int      info;
01244 #ifdef USE_FLOAT
01245     ssyevr_(&jobz,      // (IN) JOBZ:  'V'=compute eigenvals and eigenvecs
01246             &range,     // (IN) RANGE: 'A'=all eigenvals
01247             &uplo,      // (IN) UPLO:  'L'=lower triangle of A is stored
01248             &N,         // (IN) N
01249             A.GetPtr(), // (IN/OUT) float  * A
01250             &N,         // (IN) LDA
01251             &vl,        // (IN) VL for RANGE='V' ** not used **
01252             &vu,        // (IN) VU for RANGE='V' ** not used **
01253             &il,        // (IN) IL for RANGE='I' ** not used **
01254             &iu,        // (IN) IU for RANGE='I' ** not used **
01255             &abstol,    // (IN) ABSTOL
01256             &m,         // (OUT) M = total number of eigenvals found (M=N)
01257             L.GetPtr(), // (OUT) sorted (ascending) eigenvalues
01258             Q.GetPtr(), // (OUT) columns contain the eigenvecs corresp. to eigenvs
01259             &N,         // (IN) LDQ
01260             isuppz,     // (OUT) int    * ISUPPZ(2*max(1,M))
01261             work,       // (OUT) float  * WORK(LWORK)
01262             &lwork,     // (IN)  int    * LWORK >= max(1,26*N)
01263             iwork,      // (OUT) int    * IWORK(LIWORK)
01264             &liwork,    // (IN)  int    * LIWORK >= max(1,10*N)
01265             &info);     // (OUT) INFO
01266 #else
01267     dsyevr_(&jobz,      // (IN) JOBZ:  'V'=compute eigenvals and eigenvecs
01268             &range,     // (IN) RANGE: 'A'=all eigenvals
01269             &uplo,      // (IN) UPLO:  'L'=lower triangle of A is stored
01270             &N,         // (IN) N
01271             A.GetPtr(), // (IN/OUT) double * A
01272             &N,         // (IN) LDA
01273             &vl,        // (IN) VL for RANGE='V' ** not used **
01274             &vu,        // (IN) VU for RANGE='V' ** not used **
01275             &il,        // (IN) IL for RANGE='I' ** not used **
01276             &iu,        // (IN) IU for RANGE='I' ** not used **
01277             &abstol,    // (IN) ABSTOL
01278             &m,         // (OUT) M = total number of eigenvals found (M=N)
01279             L.GetPtr(), // (OUT) sorted (ascending) eigenvalues
01280             Q.GetPtr(), // (OUT) columns contain the eigenvecs corresp. to eigenvs
01281             &N,         // (IN) LDQ
01282             isuppz,     // (OUT) int    * ISUPPZ(2*max(1,M))
01283             work,       // (OUT) double * WORK(LWORK)
01284             &lwork,     // (IN)  int    * LWORK >= max(1,26*N)
01285             iwork,      // (OUT) int    * IWORK(LIWORK)
01286             &liwork,    // (IN)  int    * LIWORK >= max(1,10*N)
01287             &info);     // (OUT) INFO
01288 #endif // #ifdef USE_FLOAT
01289     delete [] iwork;
01290     delete [] work;
01291     delete [] isuppz;
01292     return info;
01293 }
01294 
01295 // }}}
01296 
01297 // {{{ Derived
01303 // ###################################################################################################################### DERIVED
01304 
01305 // -------------------------------------------------------------------------- Syep
01307 inline int Syep(Matrix<REAL> & A, 
01308                 Vector<REAL> & L, 
01309                 Matrix<REAL> P[]) 
01310 {
01311     int N = A.Rows();
01312     assert(A.Cols()==N);
01313     assert(L.Size()==N);
01314     int info;
01315     Matrix<REAL> Q(N,N); // columns are eigenVector<REAL>s
01316     info = Syevr(A,L,Q);
01317     if (info==0)
01318     {
01319         for (int i=0; i<N; ++i)
01320         {
01321             assert(P[i].Rows()==N); assert(P[i].Cols()==N);
01322             P[i].SetValues(0.0);
01323 #ifdef USE_FLOAT
01324             throw new Fatal(_("LinAlg::Syep not yet."));
01325             //cblas_sger(CblasColMajor, // CBLAS_ORDER
01326                        //N,             // M
01327                        //N,             // N
01328                        //1,             // Alpha
01329                        //&(Q(0,i)),     // const float  * X
01330                        //1,             // incX
01331                        //&(Q(0,i)),     // const float  * Y
01332                        //1,             // incY
01333                        //P[i].GetPtr(), // float  * A
01334                        //N);            // LDA
01335 #else
01336             throw new Fatal(_("LinAlg::Syep not yet."));
01337             //cblas_dger(CblasColMajor, // CBLAS_ORDER
01338                        //N,             // M
01339                        //N,             // N
01340                        //1,             // Alpha
01341                        //&(Q(0,i)),     // const double * X
01342                        //1,             // incX
01343                        //&(Q(0,i)),     // const double * Y
01344                        //1,             // incY
01345                        //P[i].GetPtr(), // double * A
01346                        //N);            // LDA
01347 #endif // #ifdef USE_FLOAT
01348         }
01349     }
01350     return info;
01351 }
01352 
01353 // -------------------------------------------------------------------------- Syif
01355 inline int Syif(Matrix<REAL> & A   , 
01356                 Matrix<REAL> & IF  , 
01357                 REAL (*func) (REAL)) 
01358 {
01359     int N = A.Rows();
01360     assert(A.Cols()==N);
01361     assert(IF.Rows()==N); assert(IF.Cols()==N);
01362     int    info;
01363     REAL   alpha;
01364     Matrix<REAL> Q(N,N); // columns are eigenVector<REAL>s
01365     Vector<REAL> L(N);   // components are eigenvalues
01366     info = Syevr(A,L,Q);
01367     if (info==0)
01368     {
01369         IF.SetValues(0.0);
01370         for (int i=0; i<N; ++i)
01371         {
01372             alpha = (*func) (L(i));
01373 #ifdef USE_FLOAT
01374             throw new Fatal(_("LinAlg::Syif not yet."));
01375             //cblas_sger(CblasColMajor, // CBLAS_ORDER
01376                        //N,             // M
01377                        //N,             // N
01378                        //alpha,         // Alpha
01379                        //&(Q(0,i)),     // const float  * X
01380                        //1,             // incX
01381                        //&(Q(0,i)),     // const float  * Y
01382                        //1,             // incY
01383                        //&(IF(0,0)),    // float  * A
01384                        //N);            // LDA
01385 #else
01386             throw new Fatal(_("LinAlg::Syif not yet."));
01387             //cblas_dger(CblasColMajor, // CBLAS_ORDER
01388                        //N,             // M
01389                        //N,             // N
01390                        //alpha,         // Alpha
01391                        //&(Q(0,i)),     // const double * X
01392                        //1,             // incX
01393                        //&(Q(0,i)),     // const double * Y
01394                        //1,             // incY
01395                        //&(IF(0,0)),    // double * A
01396                        //N);            // LDA
01397 #endif // #ifdef USE_FLOAT
01398         }
01399     }
01400     return info;
01401 }
01402  // end of group Derived
01404 // }}}
01405 
01406 #endif
01407  // end of group LAwrap
01409 
01410 }; // namespace LinAlg
01411 
01412 #endif //#define MECHSYS_LINALG_LAWRAP_H
01413 
01414 // vim:fdm=marker