MFEM v2.0
densemat.hpp
Go to the documentation of this file.
00001 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
00002 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
00003 // reserved. See file COPYRIGHT for details.
00004 //
00005 // This file is part of the MFEM library. For more information and source code
00006 // availability see http://mfem.googlecode.com.
00007 //
00008 // MFEM is free software; you can redistribute it and/or modify it under the
00009 // terms of the GNU Lesser General Public License (as published by the Free
00010 // Software Foundation) version 2.1 dated February 1999.
00011 
00012 #ifndef MFEM_DENSEMAT
00013 #define MFEM_DENSEMAT
00014 
00016 class DenseMatrix : public Matrix
00017 {
00018    friend class DenseTensor;
00019 
00020 private:
00021    int height;
00022    double *data;
00023 
00024    friend class DenseMatrixInverse;
00025    friend void Mult(const DenseMatrix &b,
00026                     const DenseMatrix &c,
00027                     DenseMatrix &a);
00028 
00029    void Eigensystem(Vector &ev, DenseMatrix *evect = NULL);
00030 
00031 public:
00034    DenseMatrix();
00035 
00037    DenseMatrix(const DenseMatrix &);
00038 
00040    explicit DenseMatrix(int s);
00041 
00043    DenseMatrix(int m, int n);
00044 
00046    DenseMatrix(const DenseMatrix &mat, char ch);
00047 
00048    DenseMatrix(double *d, int h, int w) : Matrix(w) { data = d; height = h; }
00049    void UseExternalData(double *d, int h, int w)
00050    { data = d; height = h; size = w; }
00051 
00052    void ClearExternalData() { data = NULL; height = size = 0; }
00053 
00055    void SetSize(int s);
00056 
00058    void SetSize(int h, int w);
00059 
00061    inline int Height() const { return height; }
00062 
00064    inline int Width() const { return size; }
00065 
00067    inline double *Data() const { return data; }
00068 
00070    inline double &operator()(int i, int j);
00071 
00073    inline const double &operator()(int i, int j) const;
00074 
00076    double operator*(const DenseMatrix &m) const;
00077 
00079    double Trace() const;
00080 
00082    virtual double &Elem(int i, int j);
00083 
00085    virtual const double &Elem(int i, int j) const;
00086 
00088    void Mult(const double *x, double *y) const;
00089 
00091    virtual void Mult(const Vector &x, Vector &y) const;
00092 
00094    virtual void MultTranspose(const Vector &x, Vector &y) const;
00095 
00097    void AddMult(const Vector &x, Vector &y) const;
00098 
00100    double InnerProduct(const double *x, const double *y) const;
00101 
00103    double InnerProduct(const Vector &x, const Vector &y) const
00104    { return InnerProduct((const double *)x, (const double *)y); }
00105 
00107    virtual MatrixInverse *Inverse() const;
00108 
00110    void Invert();
00111 
00113    double Det() const;
00114 
00115    double Weight() const;
00116 
00118    void Add(const double c, DenseMatrix &A);
00119 
00121    DenseMatrix &operator=(double c);
00122 
00124    DenseMatrix &operator=(const DenseMatrix &m);
00125 
00126    DenseMatrix &operator+=(DenseMatrix &m);
00127 
00128    DenseMatrix &operator-=(DenseMatrix &m);
00129 
00130    DenseMatrix &operator*=(double c);
00131 
00133    void Neg();
00134 
00136    void Norm2(double *v) const;
00137 
00138    double FNorm() const;
00139 
00140    void Eigenvalues(Vector &ev)
00141    { Eigensystem(ev); }
00142 
00143    void Eigenvalues(Vector &ev, DenseMatrix &evect)
00144    { Eigensystem(ev, &evect); }
00145 
00146    void Eigensystem(Vector &ev, DenseMatrix &evect)
00147    { Eigensystem(ev, &evect); }
00148 
00149    void SingularValues(Vector &sv) const;
00150 
00152    double CalcSingularvalue(const int i) const;
00153 
00156    void CalcEigenvalues(double *lambda, double *vec) const;
00157 
00158    void GetColumn(int c, Vector &col);
00159 
00160    void GetColumnReference(int c, Vector &col)
00161    { col.SetDataAndSize(data + c * height, height); }
00162 
00164    void Diag(double c, int n);
00166    void Diag(double *diag, int n);
00167 
00169    void Transpose();
00171    void Transpose(DenseMatrix &A);
00173    void Symmetrize();
00174 
00175    void Lump();
00176 
00181    void GradToCurl(DenseMatrix &curl);
00186    void GradToDiv(Vector &div);
00187 
00189    void CopyRows(DenseMatrix &A, int row1, int row2);
00191    void CopyCols(DenseMatrix &A, int col1, int col2);
00193    void CopyMN(DenseMatrix &A, int m, int n, int Aro, int Aco);
00195    void CopyMN(DenseMatrix &A, int row_offset, int col_offset);
00197    void CopyMNt(DenseMatrix &A, int row_offset, int col_offset);
00199    void CopyMNDiag(double c, int n, int row_offset, int col_offset);
00201    void CopyMNDiag(double *diag, int n, int row_offset, int col_offset);
00202 
00204    void AddMatrix(DenseMatrix &A, int ro, int co);
00206    void AddMatrix(double a, DenseMatrix &A, int ro, int co);
00207 
00209    void AddToVector(int offset, Vector &v) const;
00211    void GetFromVector(int offset, const Vector &v);
00214    void AdjustDofDirection(Array<int> &dofs);
00215 
00217    virtual void Print(ostream &out = cout, int width = 4) const;
00219    virtual void PrintT(ostream &out = cout, int width = 4) const;
00220 
00222    void TestInversion();
00223 
00225    virtual ~DenseMatrix();
00226 };
00227 
00229 void Add(const DenseMatrix &A, const DenseMatrix &B,
00230          double alpha, DenseMatrix &C);
00231 
00233 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
00234 
00236 void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja);
00237 
00239 void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat);
00240 
00243 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
00244 
00246 void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);
00247 
00249 void MultAAt(const DenseMatrix &a, DenseMatrix &aat);
00250 
00252 void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);
00253 
00255 void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);
00256 
00258 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
00259 
00261 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
00262 
00264 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);
00265 
00267 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
00268 
00270 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
00271 
00273 void MultVVt(const Vector &v, DenseMatrix &vvt);
00274 
00275 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
00276 
00278 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
00279 
00281 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
00282 
00283 
00286 class DenseMatrixInverse : public MatrixInverse
00287 {
00288 private:
00289    double *data;
00290 #ifdef MFEM_USE_LAPACK
00291    int *ipiv;
00292 #endif
00293 
00294 public:
00297    DenseMatrixInverse(const DenseMatrix &mat);
00298 
00300    DenseMatrixInverse(const DenseMatrix *mat);
00301 
00303    void Factor();
00304 
00306    void Factor(const DenseMatrix &mat);
00307 
00309    virtual void Mult(const Vector &x, Vector &y) const;
00310 
00312    virtual ~DenseMatrixInverse();
00313 };
00314 
00315 
00316 class DenseMatrixEigensystem
00317 {
00318    DenseMatrix &mat;
00319 
00320    Vector      EVal;
00321    DenseMatrix EVect;
00322 
00323    Vector ev;
00324 
00325    double *work;
00326    char jobz, uplo;
00327    int n, lwork, info;
00328 
00329 public:
00330 
00331    DenseMatrixEigensystem(DenseMatrix &m);
00332    void Eval();
00333    double Eigenvalue(int i) { return EVal(i); }
00334    const Vector &Eigenvector(int i)
00335    {
00336       ev.SetData(EVect.Data() + i * EVect.Height());
00337       return ev;
00338    }
00339    ~DenseMatrixEigensystem();
00340 };
00341 
00342 
00343 class DenseMatrixSVD
00344 {
00345    Vector sv;
00346 
00347    int m, n;
00348    char jobu, jobvt;
00349    int lwork, info;
00350    double *work;
00351 
00352    void Init();
00353 public:
00354 
00355    DenseMatrixSVD(DenseMatrix &M);
00356    DenseMatrixSVD(int h, int w);
00357    void Eval(DenseMatrix &M);
00358    double Singularvalue(int i) { return sv(i); }
00359    ~DenseMatrixSVD();
00360 };
00361 
00362 
00364 class DenseTensor
00365 {
00366 private:
00367    DenseMatrix Mk;
00368    double *tdata;
00369    int nk;
00370 
00371 public:
00372    DenseTensor() { nk = 0; tdata = NULL; }
00373 
00374    DenseTensor(int i, int j, int k)
00375       : Mk(NULL, i, j)
00376    { nk = k; tdata = new double[i*j*k]; }
00377 
00378    int SizeI() const { return Mk.Height(); }
00379    int SizeJ() const { return Mk.Width(); }
00380    int SizeK() const { return nk; }
00381 
00382    void SetSize(int i, int j, int k)
00383    {
00384       delete [] tdata;
00385       Mk.UseExternalData(NULL, i, j);
00386       nk = k;
00387       tdata = new double[i*j*k];
00388    }
00389 
00390    DenseMatrix &operator()(int k) { Mk.data = GetData(k); return Mk; }
00391 
00392    double *GetData(int k) { return tdata+k*Mk.Height()*Mk.Width(); }
00393 
00394    double *Data() { return tdata; }
00395 
00396    ~DenseTensor() { delete [] tdata; Mk.ClearExternalData(); }
00397 };
00398 
00399 
00400 // Inline methods
00401 
00402 inline double &DenseMatrix::operator()(int i, int j)
00403 {
00404 #ifdef MFEM_DEBUG
00405    if ( data == 0 || i < 0 || i >= height || j < 0 || j >= size )
00406       mfem_error("DenseMatrix::operator()");
00407 #endif
00408 
00409    return data[i+j*height];
00410 }
00411 
00412 inline const double &DenseMatrix::operator()(int i, int j) const
00413 {
00414 #ifdef MFEM_DEBUG
00415    if ( data == 0 || i < 0 || i >= height || j < 0 || j >= size )
00416       mfem_error("DenseMatrix::operator() const");
00417 #endif
00418 
00419    return data[i+j*height];
00420 }
00421 
00422 #endif
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines