portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-15 02:05:37 +00:00
Code Releases Activity
dc5b909f99
Grid/lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h.bak2

1042 lines
30 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Chulwoo Jung
Author: Yong-Chull Jang <ypj@quark.phy.bnl.gov>
Author: Guido Cossu
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef GRID_IRBL_H
#define GRID_IRBL_H
#include <string.h> //memset
#define clog std::cout << GridLogMessage
namespace Grid {
/////////////////////////////////////////////////////////////
// Implicitly restarted block lanczos
/////////////////////////////////////////////////////////////
template<class Field>
class ImplicitlyRestartedBlockLanczos {
private:
std::string cname = std::string("ImplicitlyRestartedBlockLanczos");
int MaxIter; // Max iterations
int Nstop; // Number of evecs checked for convergence
int Nu; // Numbeer of vecs in the unit block
int Nk; // Number of converged sought
int Nm; // total number of vectors
int Nblock_k; // Nk/Nu
int Nblock_m; // Nm/Nu
RealD eresid;
IRLdiagonalisation diagonalisation;
////////////////////////////////////
// Embedded objects
////////////////////////////////////
SortEigen<Field> _sort;
LinearOperatorBase<Field> &_Linop;
OperatorFunction<Field> &_poly;
/////////////////////////
// Constructor
/////////////////////////
public:
ImplicitlyRestartedBlockLanczos(LinearOperatorBase<Field> &Linop, // op
OperatorFunction<Field> & poly, // polynomial
int _Nstop, // really sought vecs
int _Nu, // vecs in the unit block
int _Nk, // sought vecs
int _Nm, // total vecs
RealD _eresid, // resid in lmd deficit
int _MaxIter, // Max iterations
IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen)
: _Linop(Linop), _poly(poly),
Nstop(_Nstop), Nu(_Nu), Nk(_Nk), Nm(_Nm),
Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu),
//eresid(_eresid), MaxIter(10),
eresid(_eresid), MaxIter(_MaxIter),
diagonalisation(_diagonalisation)
{ assert( (Nk%Nu==0) && (Nm%Nu==0) ); };
////////////////////////////////
// Helpers
////////////////////////////////
static RealD normalize(Field& v)
{
RealD nn = norm2(v);
nn = sqrt(nn);
v = v * (1.0/nn);
return nn;
}
void orthogonalize(Field& w, std::vector<Field>& evec, int k)
{
typedef typename Field::scalar_type MyComplex;
MyComplex ip;
for(int j=0; j<k; ++j){
ip = innerProduct(evec[j],w);
w = w - ip * evec[j];
}
normalize(w);
}
void orthogonalize_blockhead(Field& w, std::vector<Field>& evec, int k, int Nu)
{
typedef typename Field::scalar_type MyComplex;
MyComplex ip;
for(int j=0; j<k; ++j){
ip = innerProduct(evec[j*Nu],w);
w = w - ip * evec[j*Nu];
}
normalize(w);
}
/* Rudy Arthur's thesis pp.137
------------------------
Require: M > K P = M − K †
Compute the factorization AVM = VM HM + fM eM
repeat
Q=I
for i = 1,...,P do
QiRi =HM −θiI Q = QQi
H M = Q †i H M Q i
end for
βK =HM(K+1,K) σK =Q(M,K)
r=vK+1βK +rσK
VK =VM(1:M)Q(1:M,1:K)
HK =HM(1:K,1:K)
→AVK =VKHK +fKe†K † Extend to an M = K + P step factorization AVM = VMHM + fMeM
until convergence
*/
void calc(std::vector<RealD>& eval,
std::vector<Field>& evec,
const std::vector<Field>& src, int& Nconv)
{
std::string fname = std::string(cname+"::calc()");
GridBase *grid = evec[0]._grid;
assert(grid == src[0]._grid);
assert( Nu = src.size() );
clog << std::string(74,'*') << std::endl;
clog << fname + " starting iteration 0 / "<< MaxIter<< std::endl;
clog << std::string(74,'*') << std::endl;
clog <<" -- seek Nk = "<< Nk <<" vectors"<< std::endl;
clog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl;
clog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl;
clog <<" -- size of eval = "<< eval.size() << std::endl;
clog <<" -- size of evec = "<< evec.size() << std::endl;
if ( diagonalisation == IRLdiagonaliseWithEigen ) {
clog << "Diagonalisation is Eigen "<< std::endl;
} else {
abort();
}
clog << std::string(74,'*') << std::endl;
assert(Nm == evec.size() && Nm == eval.size());
std::vector<std::vector<ComplexD>> lmd(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<std::vector<ComplexD>> lme(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<std::vector<ComplexD>> lmd2(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<std::vector<ComplexD>> lme2(Nu,std::vector<ComplexD>(Nm,0.0));
std::vector<RealD> eval2(Nm);
Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
std::vector<int> Iconv(Nm);
std::vector<Field> B(Nm,grid); // waste of space replicating
std::vector<Field> f(Nu,grid);
std::vector<Field> f_copy(Nu,grid);
Field v(grid);
Nconv = 0;
RealD beta_k;
// set initial vector
for (int i=0; i<Nu; ++i) {
clog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
evec[i] = src[i];
orthogonalize(evec[i],evec,i);
clog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
}
// initial Nblock_k steps
for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
// restarting loop begins
int iter;
for(iter = 0; iter<MaxIter; ++iter){
clog <<" **********************"<< std::endl;
clog <<" Restart iteration = "<< iter << std::endl;
clog <<" **********************"<< std::endl;
// additional (Nblock_m - Nblock_k) steps
for(int b=Nblock_k; b<Nblock_m; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
//for(int k=0; k<Nm; ++k) {
// clog << "ckpt A1: lme[" << k << "] = " << lme[0][k] << '\n';
//}
//for(int k=0; k<Nm; ++k) {
// clog << "ckpt A2: lmd[" << k << "] = " << lmd[0][k] << '\n';
//}
// getting eigenvalues
for(int u=0; u<Nu; ++u){
for(int k=0; k<Nm; ++k){
lmd2[u][k] = lmd[u][k];
lme2[u][k] = lme[u][k];
}
}
Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
//for(int k=0; k<Nm; ++k){
// clog << "ckpt D " << '\n';
// clog << "eval2 [" << k << "] = " << eval2[k] << std::endl;
//}
// sorting
_sort.push(eval2,Nm);
//for(int k=0; k<Nm; ++k){
// clog << "ckpt E " << '\n';
// clog << "eval2 [" << k << "] = " << eval2[k] << std::endl;
//}
// Implicitly shifted QR transformations
Eigen::MatrixXcd BTDM = Eigen::MatrixXcd::Identity(Nm,Nm);
Q = Eigen::MatrixXcd::Identity(Nm,Nm);
unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
for(int ip=Nk; ip<Nm; ++ip){
//clog << "ckpt B1: shift[" << ip << "] = " << eval2[ip] << endl;
//for (int i=0; i<Nm; ++i) {
// for (int j=0; j<Nm; ++j) {
// clog << "ckpt B2: M(" << ip << ")[" << i << "," << j << "] = " << BTDM(i,j) << '\n';
// }
//}
//for (int i=0; i<Nm; ++i) {
// for (int j=0; j<Nm; ++j) {
// clog << "ckpt B3: Q(" << ip << ")[" << i << "," << j << "] = " << Q(i,j) << '\n';
// }
//}
shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q);
}
//BTDM = Q.adjoint()*(BTDM*Q);
//for (int i=0; i<Nm; ++i ) {
// for (int j=i+1; j<Nm; ++j ) {
// BTDM(i,j) = conj(BTDM(j,i));
// }
// //BTDM(i,i) = real(BTDM(i,i));
//}
//for (int i=0; i<Nm; ++i) {
// for (int j=0; j<Nm; ++j) {
// clog << "ckpt B2: M(" << Nm-Nk << ")[" << i << "," << j << "] = " << BTDM(i,j) << '\n';
// }
//}
//for (int i=0; i<Nm; ++i) {
// for (int j=0; j<Nm; ++j) {
// clog << "ckpt B3: Q(" << Nm-Nk << ")[" << i << "," << j << "] = " << Q(i,j) << '\n';
// }
//}
packHermitBlockTriDiagMatfromEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
//for (int i=0; i<Nm; ++i) {
// clog << "ckpt C1: lme[" << i << "] = " << lme[0][i] << '\n';
//}
//for (int i=0; i<Nm; ++i) {
// clog << "ckpt C2: lmd[" << i << "] = " << lmd[0][i] << '\n';
//}
for(int i=0; i<Nk+Nu; ++i) B[i] = 0.0;
for(int j=0; j<Nk+Nu; ++j){
for(int k=0; k<Nm; ++k){
B[j].checkerboard = evec[k].checkerboard;
B[j] += evec[k]*Q(k,j);
}
}
for(int i=0; i<Nk+Nu; ++i) {
evec[i] = B[i];
//clog << "ckpt F: norm2_evec[= " << i << "]" << norm2(evec[i]) << std::endl;
}
// residual vector
#if 0 // ypj[fixme] temporary to check a case when block has one vector
for ( int i=0; i<Nu; ++i) f_copy[i] = f[i];
for ( int i=0; i<Nu; ++i) {
f[i] = f_copy[0]*lme[0][Nm-Nu+i];
for ( int j=1; j<Nu; ++j) {
f[i] += f_copy[j]*lme[j][Nm-Nu+i];
}
//clog << "ckpt C (i= " << i << ")" << '\n';
//clog << "norm2(f) = " << norm2(f[i]) << std::endl;
}
// ypj[fixme] temporary to check a case when block has one vector
// Compressed vector f and beta(k2)
f[0] *= Q(Nm-1,Nk-1);
f[0] += lme[0][Nk-1] * evec[Nk]; // was commented out
std::cout<< GridLogMessage<<"ckpt D1: Q[Nm-1,Nk-1] = "<<Q(Nm-1,Nk-1)<<std::endl;
beta_k = norm2(f[0]);
beta_k = sqrt(beta_k);
std::cout<< GridLogMessage<<"ckpt D2: beta(k) = "<<beta_k<<std::endl;
RealD betar = 1.0/beta_k;
evec[Nk] = betar * f[0];
lme[0][Nk-1] = beta_k;
#else
blockwiseStep(lmd,lme,evec,f,f_copy,Nblock_k-1);
// set initial vector
//for (int i=0; i<Nu; ++i) {
// orthogonalize(evec[i],evec,i);
//}
//// initial Nblock_k steps
//for(int b=0; b<Nblock_k; ++b) blockwiseStep(lmd,lme,evec,f,f_copy,b);
#endif
// Convergence test
for(int u=0; u<Nu; ++u){
for(int k=0; k<Nm; ++k){
lmd2[u][k] = lmd[u][k];
lme2[u][k] = lme[u][k];
}
}
Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
for(int k = 0; k<Nk; ++k) B[k]=0.0;
for(int j = 0; j<Nk; ++j){
for(int k = 0; k<Nk; ++k){
B[j].checkerboard = evec[k].checkerboard;
B[j] += evec[k]*Qt(k,j);
}
}
//for (int k=0; k<Nk; ++k) {
// orthogonalize(B[k],B,k);
//}
//for (int i=0; i<Nk; ++i) {
// for (int j=0; j<Nk; ++j) {
// clog << "ckpt H1: R[" << i << "," << j << "] = " << Qt(i,j) << '\n';
// }
//}
//for (int i=0; i<Nk; ++i) {
// clog << "ckpt H2: eval2[" << i << "] = " << eval2[i] << '\n';
//}
//for(int j=0; j<Nk; ++j) {
// clog << "ckpt I: norm2_B[ " << j << "]" << norm2(B[j]) << std::endl;
//}
Nconv = 0;
for(int i=0; i<Nk; ++i){
_Linop.HermOp(B[i],v);
RealD vnum = real(innerProduct(B[i],v)); // HermOp.
RealD vden = norm2(B[i]);
eval2[i] = vnum/vden;
v -= eval2[i]*B[i];
RealD vv = norm2(v);
std::cout.precision(13);
clog << "[" << std::setw(3)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
std::cout << "eval = "<<std::setw(25)<< std::setiosflags(std::ios_base::left)<< eval2[i];
std::cout << " |H B[i] - eval[i]B[i]|^2 "<< std::setw(25)<< std::setiosflags(std::ios_base::right)<< vv<< std::endl;
// change the criteria as evals are supposed to be sorted, all evals smaller(larger) than Nstop should have converged
if( (vv<eresid*eresid) && (i == Nconv) ){
//if( (vv<eresid*eresid) ){
Iconv[Nconv] = i;
++Nconv;
}
} // i-loop end
clog <<" #modes converged: "<<Nconv<<std::endl;
if( Nconv>=Nstop ){
goto converged;
}
} // end of iter loop
clog <<"**************************************************************************"<< std::endl;
std::cout<< GridLogError << fname + " NOT converged.";
clog <<"**************************************************************************"<< std::endl;
abort();
converged:
// Sorting
eval.resize(Nconv);
evec.resize(Nconv,grid);
for(int i=0; i<Nconv; ++i){
eval[i] = eval2[Iconv[i]];
evec[i] = B[Iconv[i]];
}
_sort.push(eval,evec,Nconv);
clog <<"**************************************************************************"<< std::endl;
clog << fname + " CONVERGED ; Summary :\n";
clog <<"**************************************************************************"<< std::endl;
clog << " -- Iterations = "<< iter << "\n";
//clog << " -- beta(k) = "<< beta_k << "\n";
clog << " -- Nconv = "<< Nconv << "\n";
clog <<"**************************************************************************"<< std::endl;
}
private:
/* Saad PP. 195
1. Choose an initial vector v1 of 2-norm unity. Set β1 ≡ 0, v0 ≡ 0
2. For k = 1,2,...,m Do:
3. wk:=Avk−βkv_{k−1}
4. αk:=(wk,vk) //
5. wk:=wk−αkvk // wk orthog vk
6. βk+1 := ∥wk∥2. If βk+1 = 0 then Stop
7. vk+1 := wk/βk+1
8. EndDo
*/
void blockwiseStep(std::vector<std::vector<ComplexD>>& lmd,
std::vector<std::vector<ComplexD>>& lme,
std::vector<Field>& evec,
std::vector<Field>& w,
std::vector<Field>& w_copy,
int b)
{
const RealD tiny = 1.0e-20;
int Nu = w.size();
int Nm = evec.size();
assert( b < Nm/Nu );
// converts block index to full indicies for an interval [L,R)
int L = Nu*b;
int R = Nu*(b+1);
Real beta;
// 3. wk:=Avk−βkv_{k−1}
for (int k=L, u=0; k<R; ++k, ++u) {
_poly(_Linop,evec[k],w[u]);
}
if (b>0) {
for (int u=0; u<Nu; ++u) {
//for (int k=L-Nu; k<L; ++k) {
for (int k=L-Nu+u; k<L; ++k) {
w[u] = w[u] - evec[k] * conjugate(lme[u][k]);
//clog << "ckpt A (k= " << k+1 << ")" << '\n';
//clog << "lme = " << lme[u][k] << '\n';
//clog << "lme = " << conjugate(lme[u][k]) << '\n';
}
//clog << "norm(w) = " << norm2(w[u]) << std::endl;
}
}
// 4. αk:=(vk,wk)
//for (int u=0; u<Nu; ++u) {
// for (int k=L; k<R; ++k) {
// lmd[u][k] = innerProduct(evec[k],w[u]); // lmd = transpose of alpha
// }
// lmd[u][L+u] = real(lmd[u][L+u]); // force diagonal to be real
//}
for (int u=0; u<Nu; ++u) {
for (int k=L+u; k<R; ++k) {
lmd[u][k] = innerProduct(evec[k],w[u]); // lmd = transpose of alpha
lmd[k-L][u+L] = conjugate(lmd[u][k]); // force hermicity
}
lmd[u][L+u] = real(lmd[u][L+u]); // force diagonal to be real
}
//for (int u=0; u<Nu; ++u) {
// for (int k=L; k<L+u+1; ++k) {
// clog << "ckpt B: alphaT[" << u+L << "," << k << "] = " << lmd[u][k]
// << ", alphaT - conj(alphaT) = " << lmd[u][k] - conjugate(lmd[k-L][u+L]) << std::endl;
// }
//}
// 5. wk:=wk−αkvk
for (int u=0; u<Nu; ++u) {
for (int k=L; k<R; ++k) {
w[u] = w[u] - evec[k]*lmd[u][k];
}
w_copy[u] = w[u];
}
// In block version, the steps 6 and 7 in Lanczos construction is
// replaced by the QR decomposition of new basis block.
// It results block version beta and orthonormal block basis.
// Here, QR decomposition is done by using Gram-Schmidt.
for (int u=0; u<Nu; ++u) {
for (int k=L; k<R; ++k) {
lme[u][k] = 0.0;
}
}
#if 0
beta = normalize(w[0]);
for (int u=1; u<Nu; ++u) {
//orthogonalize(w[u],w_copy,u);
orthogonalize(w[u],w,u);
}
#else
// re-orthogonalization for numerical stability
for (int u=0; u<Nu; ++u) {
orthogonalize(w[u],evec,R);
}
// QR part
for (int u=1; u<Nu; ++u) {
orthogonalize(w[u],w,u);
}
#endif
for (int u=0; u<Nu; ++u) {
//for (int v=0; v<Nu; ++v) {
for (int v=u; v<Nu; ++v) {
lme[u][L+v] = innerProduct(w[u],w_copy[v]);
}
lme[u][L+u] = real(lme[u][L+u]); // force diagonal to be real
}
//lme[0][L] = beta;
for (int u=0; u<Nu; ++u) {
clog << "norm2(w[" << u << "])= "<< norm2(w[u]) << std::endl;
for (int k=L+u; k<R; ++k) {
clog <<" In block "<< b << ",";
std::cout <<" beta[" << u << "," << k-L << "] = ";
std::cout << lme[u][k] << std::endl;
}
}
#if 0
// re-orthogonalization for numerical stability
if (b>0) {
for (int u=0; u<Nu; ++u) {
orthogonalize(w[u],evec,R);
}
for (int u=1; u<Nu; ++u) {
orthogonalize(w[u],w,u);
}
}
//if (b>0) {
// orthogonalize_blockhead(w[0],evec,b,Nu);
// for (int u=1; u<Nu; ++u) {
// orthogonalize(w[u],w,u);
// }
//}
#endif
if (b < Nm/Nu-1) {
for (int u=0; u<Nu; ++u) {
evec[R+u] = w[u];
}
}
}
void diagonalize_Eigen(std::vector<RealD>& eval,
std::vector<std::vector<ComplexD>>& lmd,
std::vector<std::vector<ComplexD>>& lme,
int Nu, int Nk, int Nm,
Eigen::MatrixXcd & Qt, // Nm x Nm
GridBase *grid)
{
assert( Nk%Nu == 0 && Nm%Nu == 0 );
assert( Nk <= Nm );
Eigen::MatrixXcd BlockTriDiag = Eigen::MatrixXcd::Zero(Nk,Nk);
for ( int u=0; u<Nu; ++u ) {
for (int k=0; k<Nk; ++k ) {
BlockTriDiag(k,u+(k/Nu)*Nu) = lmd[u][k];
}
}
for ( int u=0; u<Nu; ++u ) {
for (int k=Nu; k<Nk; ++k ) {
BlockTriDiag(k-Nu,u+(k/Nu)*Nu) = conjugate(lme[u][k-Nu]);
BlockTriDiag(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
}
}
//std::cout << BlockTriDiag << std::endl;
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver(BlockTriDiag);
for (int i = 0; i < Nk; i++) {
eval[Nk-1-i] = eigensolver.eigenvalues()(i);
}
for (int i = 0; i < Nk; i++) {
for (int j = 0; j < Nk; j++) {
Qt(j,Nk-1-i) = eigensolver.eigenvectors()(j,i);
//Qt(Nk-1-i,j) = eigensolver.eigenvectors()(i,j);
//Qt(i,j) = eigensolver.eigenvectors()(i,j);
}
}
}
void diagonalize(std::vector<RealD>& eval,
std::vector<std::vector<ComplexD>>& lmd,
std::vector<std::vector<ComplexD>>& lme,
int Nu, int Nk, int Nm,
Eigen::MatrixXcd & Qt,
GridBase *grid)
{
Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
if ( diagonalisation == IRLdiagonaliseWithEigen ) {
diagonalize_Eigen(eval,lmd,lme,Nu,Nk,Nm,Qt,grid);
} else {
assert(0);
}
}
void unpackHermitBlockTriDiagMatToEigen(
std::vector<std::vector<ComplexD>>& lmd,
std::vector<std::vector<ComplexD>>& lme,
int Nu, int Nb, int Nk, int Nm,
Eigen::MatrixXcd& M)
{
//clog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n';
assert( Nk%Nu == 0 && Nm%Nu == 0 );
assert( Nk <= Nm );
M = Eigen::MatrixXcd::Zero(Nk,Nk);
// rearrange
for ( int u=0; u<Nu; ++u ) {
for (int k=0; k<Nk; ++k ) {
M(k,u+(k/Nu)*Nu) = lmd[u][k];
}
}
for ( int u=0; u<Nu; ++u ) {
for (int k=Nu; k<Nk; ++k ) {
M(k-Nu,u+(k/Nu)*Nu) = conjugate(lme[u][k-Nu]);
M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
}
}
//clog << "unpackHermitBlockTriDiagMatToEigen() end" << endl;
}
void packHermitBlockTriDiagMatfromEigen(
std::vector<std::vector<ComplexD>>& lmd,
std::vector<std::vector<ComplexD>>& lme,
int Nu, int Nb, int Nk, int Nm,
Eigen::MatrixXcd& M)
{
//clog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n';
assert( Nk%Nu == 0 && Nm%Nu == 0 );
assert( Nk <= Nm );
// rearrange
for ( int u=0; u<Nu; ++u ) {
for (int k=0; k<Nk; ++k ) {
lmd[u][k] = M(k,u+(k/Nu)*Nu);
}
}
for ( int u=0; u<Nu; ++u ) {
for (int k=Nu; k<Nk; ++k ) {
lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
}
}
//clog << "packHermitBlockTriDiagMatfromEigen() end" << endl;
}
#if 0
void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
RealD Dsh,
Eigen::MatrixXcd& Qprod)
{
//clog << "shiftedQRDecompEigen() begin" << '\n';
Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
Mtmp = M;
for (int i=0; i<Nm; ++i ) {
Mtmp(i,i) = M(i,i) - Dsh;
}
Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
Q = QRD.householderQ();
for (int j=0; j<Nm-2*Nu; ++j ) {
for (int i=2*Nu+j; i<Nm; ++i ) {
Q(i,j) = 0.;
}
}
M = Q.adjoint()*(M*Q);
for (int i=0; i<Nm; ++i ) {
for (int j=i+Nu; j<Nm; ++j ) {
M(i,j) = conj(M(j,i));
}
}
for (int i=0; i<Nm; ++i ) {
M(i,i) = real(M(i,i));
}
Qprod *= Q;
//clog << "shiftedQRDecompEigen() end" << endl;
}
#endif
#if 1
// assume the input matrix M is a band matrix
void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
RealD Dsh,
Eigen::MatrixXcd& Qprod)
{
//clog << "shiftedQRDecompEigen() begin" << '\n';
Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd R = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
Mtmp = M;
for (int i=0; i<Nm; ++i ) {
Mtmp(i,i) = M(i,i) - Dsh;
}
Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
Q = QRD.householderQ();
R = QRD.matrixQR(); // upper triangular part is the R matrix.
// lower triangular part used to represent series
// of Q sequence.
// equivalent operation of Qprod *= Q
//M = Eigen::MatrixXcd::Zero(Nm,Nm);
//for (int i=0; i<Nm; ++i) {
// for (int j=0; j<Nm-2*(Nu+1); ++j) {
// for (int k=0; k<2*(Nu+1)+j; ++k) {
// M(i,j) += Qprod(i,k)*Q(k,j);
// }
// }
//}
//for (int i=0; i<Nm; ++i) {
// for (int j=Nm-2*(Nu+1); j<Nm; ++j) {
// for (int k=0; k<Nm; ++k) {
// M(i,j) += Qprod(i,k)*Q(k,j);
// }
// }
//}
Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
for (int i=0; i<Nm; ++i) {
for (int j=0; j<Nm-(Nu+1); ++j) {
for (int k=0; k<Nu+1+j; ++k) {
//for (int j=0; j<Nm-2*Nu; ++j) {
// for (int k=0; k<j+2*Nu; ++k) {
Mtmp(i,j) += Qprod(i,k)*Q(k,j);
}
}
}
for (int i=0; i<Nm; ++i) {
for (int j=Nm-(Nu+1); j<Nm; ++j) {
//for (int j=Nm-2*Nu; j<Nm; ++j) {
for (int k=0; k<Nm; ++k) {
Mtmp(i,j) += Qprod(i,k)*Q(k,j);
}
}
}
//static int ntimes = 2;
//for (int j=0; j<Nm-(ntimes*Nu); ++j) {
// for (int i=ntimes*Nu+j; i<Nm; ++i) {
// Mtmp(i,j) = 0.0;
// }
//}
//ntimes++;
Qprod = Mtmp;
// equivalent operation of M = Q.adjoint()*(M*Q)
//M = Eigen::MatrixXcd::Zero(Nm,Nm);
//
//for (int a=0, i=0, kmax=0; a<Nu+1; ++a) {
// for (int j=0; j<Nm-a; ++j) {
// i = j+a;
// kmax = (Nu+1)+j;
// if (kmax > Nm) kmax = Nm;
// for (int k=i; k<kmax; ++k) {
// M(i,j) += R(i,k)*Q(k,j);
// }
// M(j,i) = conj(M(i,j));
// }
//}
//for (int i=0; i<Nm; ++i) {
// M(i,i) = real(M(i,i));
//}
M = Q.adjoint()*(M*Q);
for (int i=0; i<Nm; ++i) {
for (int j=0; j<Nm; ++j) {
if (i==j) M(i,i) = real(M(i,i));
if (j>i) M(i,j) = conj(M(j,i));
if (i-j > Nu || j-i > Nu) M(i,j) = 0.;
}
}
//clog << "shiftedQRDecompEigen() end" << endl;
}
#endif
#if 0
void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nu, int Nm,
RealD Dsh,
Eigen::MatrixXcd& Qprod)
{
//clog << "shiftedQRDecompEigen() begin" << '\n';
Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
Mtmp = M;
for (int i=0; i<Nm; ++i ) {
Mtmp(i,i) = M(i,i) - Dsh;
}
Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
Q = QRD.householderQ();
M = Q.adjoint()*(M*Q);
for (int i=0; i<Nm; ++i ) {
for (int j=i+1; j<Nm; ++j ) {
M(i,j) = conj(M(j,i));
}
}
for (int i=0; i<Nm; ++i ) {
M(i,i) = real(M(i,i));
}
#if 1
Qprod *= Q;
#else
Mtmp = Qprod*Q;
Eigen::HouseholderQR<Eigen::MatrixXcd> QRD2(Mtmp);
Qprod = QRD2.householderQ();
Mtmp -= Qprod;
clog << "Frobenius norm ||Qprod(after) - Qprod|| = " << Mtmp.norm() << std::endl;
#endif
//clog << "shiftedQRDecompEigen() end" << endl;
}
#endif
#if 0
void shiftedQRDecompEigen(Eigen::MatrixXcd& M, int Nm,
RealD Dsh,
Eigen::MatrixXcd& Qprod)
{
//clog << "shiftedQRDecompEigen() begin" << '\n';
Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
//Eigen::MatrixXcd Qtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
Mtmp = Qprod.adjoint()*(M*Qprod);
for (int i=0; i<Nm; ++i ) {
for (int j=i+1; j<Nm; ++j ) {
Mtmp(i,j) = Mtmp(j,i);
}
}
for (int i=0; i<Nm; ++i ) {
Mtmp(i,i) -= Dsh;
//Mtmp(i,i) = real(Mtmp(i,i)-Dsh);
}
Eigen::HouseholderQR<Eigen::MatrixXcd> QRD(Mtmp);
//Qtmp = Qprod*QRD.householderQ();
//Eigen::HouseholderQR<Eigen::MatrixXcd> QRD2(Qtmp);
//Qprod = QRD2.householderQ();
Qprod *= QRD.householderQ();
//ComplexD p;
//RealD r;
//r = 0.;
//for (int k=0; k<Nm; ++k) r += real(conj(Qprod(k,0))*Qprod(k,0));
//r = sqrt(r);
//for (int k=0; k<Nm; ++k) Qprod(k,0) /= r;
//
//for (int i=1; i<Nm; ++i) {
// for (int j=0; j<i; ++j) {
// p = 0.;
// for (int k=0; k<Nm; ++k) {
// p += conj(Qprod(k,j))*Qprod(k,i);
// }
// for (int k=0; k<Nm; ++k) {
// Qprod(k,i) -= p*Qprod(k,j);
// }
// }
// r = 0.;
// for (int k=0; k<Nm; ++k) r += real(conj(Qprod(k,i))*Qprod(k,i));
// r = sqrt(r);
// for (int k=0; k<Nm; ++k) Qprod(k,i) /= r;
//}
//clog << "shiftedQRDecompEigen() end" << endl;
}
#endif
void exampleQRDecompEigen(void)
{
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3,3);
Eigen::MatrixXd Q = Eigen::MatrixXd::Zero(3,3);
Eigen::MatrixXd R = Eigen::MatrixXd::Zero(3,3);
Eigen::MatrixXd P = Eigen::MatrixXd::Zero(3,3);
A(0,0) = 12.0;
A(0,1) = -51.0;
A(0,2) = 4.0;
A(1,0) = 6.0;
A(1,1) = 167.0;
A(1,2) = -68.0;
A(2,0) = -4.0;
A(2,1) = 24.0;
A(2,2) = -41.0;
clog << "matrix A before ColPivHouseholder" << std::endl;
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
}
}
clog << std::endl;
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> QRD(A);
clog << "matrix A after ColPivHouseholder" << std::endl;
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
}
}
clog << std::endl;
clog << "HouseholderQ with sequence lenth = nonzeroPiviots" << std::endl;
Q = QRD.householderQ().setLength(QRD.nonzeroPivots());
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
}
}
clog << std::endl;
clog << "HouseholderQ with sequence lenth = 1" << std::endl;
Q = QRD.householderQ().setLength(1);
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
}
}
clog << std::endl;
clog << "HouseholderQ with sequence lenth = 2" << std::endl;
Q = QRD.householderQ().setLength(2);
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
}
}
clog << std::endl;
clog << "matrixR" << std::endl;
R = QRD.matrixR();
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "R[" << i << "," << j << "] = " << R(i,j) << '\n';
}
}
clog << std::endl;
clog << "rank = " << QRD.rank() << std::endl;
clog << "threshold = " << QRD.threshold() << std::endl;
clog << "matrixP" << std::endl;
P = QRD.colsPermutation();
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "P[" << i << "," << j << "] = " << P(i,j) << '\n';
}
}
clog << std::endl;
clog << "QR decomposition without column pivoting" << std::endl;
A(0,0) = 12.0;
A(0,1) = -51.0;
A(0,2) = 4.0;
A(1,0) = 6.0;
A(1,1) = 167.0;
A(1,2) = -68.0;
A(2,0) = -4.0;
A(2,1) = 24.0;
A(2,2) = -41.0;
clog << "matrix A before Householder" << std::endl;
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
}
}
clog << std::endl;
Eigen::HouseholderQR<Eigen::MatrixXd> QRDplain(A);
clog << "HouseholderQ" << std::endl;
Q = QRDplain.householderQ();
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n';
}
}
clog << std::endl;
clog << "matrix A after Householder" << std::endl;
for ( int i=0; i<3; i++ ) {
for ( int j=0; j<3; j++ ) {
clog << "A[" << i << "," << j << "] = " << A(i,j) << '\n';
}
}
clog << std::endl;
}
};
}
#undef clog
#endif
View Git Blame Copy Permalink