PEXSI
 All Classes Namespaces Files Functions Variables Typedefs Pages
NumMat_impl.hpp
Go to the documentation of this file.
1 /*
2  Copyright (c) 2012 The Regents of the University of California,
3  through Lawrence Berkeley National Laboratory.
4 
5 Authors: Lexing Ying, Mathias Jacquelin and Lin Lin
6 
7 This file is part of PEXSI. All rights reserved.
8 
9 Redistribution and use in source and binary forms, with or without
10 modification, are permitted provided that the following conditions are met:
11 
12 (1) Redistributions of source code must retain the above copyright notice, this
13 list of conditions and the following disclaimer.
14 (2) Redistributions in binary form must reproduce the above copyright notice,
15 this list of conditions and the following disclaimer in the documentation
16 and/or other materials provided with the distribution.
17 (3) Neither the name of the University of California, Lawrence Berkeley
18 National Laboratory, U.S. Dept. of Energy nor the names of its contributors may
19 be used to endorse or promote products derived from this software without
20 specific prior written permission.
21 
22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
23 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
24 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
25 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
26 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
27 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
28 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
29 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 
33 You are under no obligation whatsoever to provide any bug fixes, patches, or
34 upgrades to the features, functionality or performance of the source code
35 ("Enhancements") to anyone; however, if you choose to make your Enhancements
36 available either publicly, or directly to Lawrence Berkeley National
37 Laboratory, without imposing a separate written license agreement for such
38 Enhancements, then you hereby grant the following license: a non-exclusive,
39 royalty-free perpetual license to install, use, modify, prepare derivative
40 works, incorporate into other computer software, distribute, and sublicense
41 such enhancements or derivative works thereof, in binary and source code form.
42  */
46 #ifndef _PEXSI_NUMMAT_IMPL_HPP_
47 #define _PEXSI_NUMMAT_IMPL_HPP_
48 
49 
50 namespace PEXSI{
51 
52  template <class F> inline void NumMat<F>::allocate(F* data) {
53  if(owndata_) {
54  if(m_>0 && n_>0) { data_ = new F[m_*n_]; if( data_ == NULL ) throw std::runtime_error("Cannot allocate memory."); } else data_=NULL;
55  if(data!=NULL){std::copy(data,data+m_*n_,data_);}
56  } else {
57  data_ = data;
58  }
59  bufsize_=m_*n_;
60  }
61  template <class F> inline void NumMat<F>::deallocate(){
62  if(owndata_) {
63  if(bufsize_>0) { delete[] data_; data_ = NULL; }
64  }
65  }
66 
67  template <class F> NumMat<F>::NumMat(Int m, Int n): m_(m), n_(n), owndata_(true) {
68  this->allocate();
69  }
70 
71  template <class F> NumMat<F>::NumMat(Int m, Int n, bool owndata, F* data): m_(m), n_(n), owndata_(owndata) {
72  this->allocate(data);
73  }
74 
75  template <class F> NumMat<F>::NumMat(const NumMat& C): m_(C.m_), n_(C.n_), owndata_(C.owndata_) {
76  this->allocate(C.data_);
77  }
78 
79  template <class F> NumMat<F>::~NumMat() {
80  this->deallocate();
81  }
82 
83  template <class F> NumMat<F>& NumMat<F>::Copy(const NumMat<F>& C) {
84  this->deallocate();
85  m_ = C.m_; n_=C.n_; owndata_=C.owndata_;
86  this->allocate(C.data_);
87  return *this;
88  }
89 
90  template <class F> NumMat<F>& NumMat<F>::operator=(const NumMat<F>& C) {
91  this->deallocate();
92  m_ = C.m_; n_=C.n_; owndata_=C.owndata_;
93  this->allocate(C.data_);
94  return *this;
95  }
96 
97  template <class F> void NumMat<F>::Resize(Int m, Int n) {
98  if( owndata_ == false ){
99  throw std::logic_error("Matrix being resized must own data.");
100  }
101 
102  if(m*n > bufsize_) {
103  this->deallocate();
104  m_ = m; n_ = n;
105  this->allocate();
106  }
107  else{
108  m_ = m; n_ = n;
109  }
110  }
111 
112  template <class F> const F& NumMat<F>::operator()(Int i, Int j) const {
113  if( i < 0 || i >= m_ ||
114  j < 0 || j >= n_ ) {
115  throw std::logic_error( "Index is out of bound." );
116  }
117  return data_[i+j*m_];
118  }
119 
120  template <class F> F& NumMat<F>::operator()(Int i, Int j) {
121  if( i < 0 || i >= m_ ||
122  j < 0 || j >= n_ ) {
123  throw std::logic_error( "Index is out of bound." );
124  }
125  return data_[i+j*m_];
126  }
127 
128  template <class F> F* NumMat<F>::VecData(Int j) const
129  {
130  if( j < 0 || j >= n_ ) {
131  throw std::logic_error( "Index is out of bound." );
132  }
133  return &(data_[j*m_]);
134  }
135 
136 
137  template <class F> inline void SetValue(NumMat<F>& M, F val)
138  {
139  std::fill(M.Data(),M.Data()+M.m()*M.n(),val);
140  }
141 
142  template <class F> inline Real Energy(const NumMat<F>& M)
143  {
144  Real sum = 0;
145  F *ptr = M.Data();
146  for (Int i=0; i < M.m()*M.n(); i++)
147  sum += abs(ptr[i]) * abs(ptr[i]);
148  return sum;
149  }
150 
151 
152  template <class F> inline void
153  Transpose ( const NumMat<F>& A, NumMat<F>& B )
154  {
155 #ifndef _RELEASE_
156  PushCallStack("Transpose");
157 #endif
158  if( A.m() != B.n() || A.n() != B.m() ){
159  B.Resize( A.n(), A.m() );
160  }
161 
162  F* Adata = A.Data();
163  F* Bdata = B.Data();
164  Int m = A.m(), n = A.n();
165 
166  for( Int i = 0; i < m; i++ ){
167  for( Int j = 0; j < n; j++ ){
168  Bdata[ j + n*i ] = Adata[ i + j*m ];
169  }
170  }
171 
172 #ifndef _RELEASE_
173  PopCallStack();
174 #endif
175 
176  return ;
177  } // ----- end of function Transpose -----
178 
179  template <class F> inline void
180  Symmetrize( NumMat<F>& A )
181  {
182 #ifndef _RELEASE_
183  PushCallStack("Symmetrize");
184 #endif
185  if( A.m() != A.n() ){
186  throw std::logic_error( "The matrix to be symmetrized should be a square matrix." );
187  }
188 
189  NumMat<F> B;
190  Transpose( A, B );
191 
192  F* Adata = A.Data();
193  F* Bdata = B.Data();
194 
195  F half = (F) 0.5;
196 
197  for( Int i = 0; i < A.m() * A.n(); i++ ){
198  *Adata = half * (*Adata + *Bdata);
199  Adata++; Bdata++;
200  }
201 
202 #ifndef _RELEASE_
203  PopCallStack();
204 #endif
205 
206  return ;
207  } // ----- end of function Symmetrize -----
208 
209 
210 } // namespace PEXSI
211 
212 #endif // _PEXSI_NUMMAT_IMPL_HPP_
Real Energy(const NumMat< F > &M)
Energy computes the L2 norm of a matrix (treated as a vector).
Definition: NumMat_impl.hpp:142
void deallocate()
Helper function freeing memory pointed by the data_ attribute.
Definition: NumMat_impl.hpp:61
void allocate(F *data=NULL)
Helper function allocating the memory pointed by the data_ attribute.
Definition: NumMat_impl.hpp:52
void SetValue(NumMat< F > &M, F val)
SetValue sets a numerical matrix to a constant val.
Definition: NumMat_impl.hpp:137
Numerical matrix.
Definition: NumMat.hpp:61