{

LAPACK_INT info = 0;

dsyev_wrapper('V','U',N,Udata,ddata,info);

return info;

{

LAPACK_INT info = 0;

std::vector<LAPACK_REAL> tau(N);

dgeqrf_wrapper(&M, &N, Qdata, &M, tau.data(), &info);

for(int i = 0; i < Rrows; i++)

for(int j = i; j < N; j++)

{

Rdata[i + j*Rrows] = Qdata[i+j*M];

}

int min = M < N ? M : N;

dorgqr_wrapper(&M, &Rrows, &min, Qdata, &M, tau.data(), &info);

return info;

{

LAPACK_INT info = 0;

std::vector<LAPACK_COMPLEX> tau(N);

zgeqrf_wrapper(&M, &N, Qdata, &M, tau.data(), &info);

for(int i = 0; i < Rrows; i++)

for(int j = i; j < N; j++)

{

Rdata[i + j*Rrows] = Qdata[i+j*M];

}

int min = M < N ? M : N;

zungqr_wrapper(&M, &Rrows, &min, Qdata, &M, tau.data(), &info);

return info;

{

LAPACK_INT info = 0;

char S = 'S';

zgesdd_wrapper(&S, &M, &N, Adata, Ddata, Udata, Vdata, &info);

return info;

{

LAPACK_INT info = 0;

char S = 'S';

dgesdd_wrapper(&S, &M, &N, Adata, Ddata, Udata, Vdata, &info);

return info;

{

LAPACK_INT info = 0;

char S = 'S';

zgesvd_wrapper(&S, &M, &N, Adata, Ddata, Udata, Vdata, &info);

return info;

{

LAPACK_INT info = 0;

char S = 'S';

dgesvd_wrapper(&S, &M, &N, Adata, Ddata, Udata, Vdata, &info);

return info;

MatrixRef const& Rr,

MatrixRef const& Ri,

MatrixRef const& Lr,

MatrixRef const& Li,

VectorRef const& dr,

VectorRef const& di)

{

auto N = ncols(M);

if(N < 1) throw std::runtime_error("diagGeneral: 0 dimensional matrix");

if(N != nrows(M))

{

printfln("M is %dx%d",nrows(M),ncols(M));

throw std::runtime_error("diagGeneral: Input Matrix must be square");

}

#ifdef DEBUG

if(!isContiguous(Rr))

throw std::runtime_error("diagGeneral: Rr must be contiguous");

if(!isContiguous(Ri))

throw std::runtime_error("diagGeneral: Ri must be contiguous");

if(Lr && !isContiguous(Lr))

throw std::runtime_error("diagGeneral: Lr must be contiguous");

if(Li && !isContiguous(Li))

throw std::runtime_error("diagGeneral: Li must be contiguous");

if(!isContiguous(dr))

throw std::runtime_error("diagGeneral: dr must be contiguous");

if(!isContiguous(di))

throw std::runtime_error("diagGeneral: di must be contiguous");

#endif

struct Diag

{

LAPACK_INT static

call(LAPACK_INT N, Real const* Mdata, Real *Ldata, Real *Rdata, Real *drdata, Real *didata)

{

auto cl = (Ldata==nullptr) ? 'N' : 'V';

return dgeev_wrapper(cl,'V',N,Mdata,drdata,didata,Ldata,Rdata);

}

LAPACK_INT static

call(LAPACK_INT N, Cplx const* Mdata, Cplx *Ldata, Cplx *Rdata, Real *drdata, Real *didata)

{

auto d = std::vector<Cplx>(N);

auto cl = (Ldata==nullptr) ? 'N' : 'V';

auto info = zgeev_wrapper(cl,'V',N,Mdata,d.data(),Ldata,Rdata);

for(size_t n = 0ul; n < d.size(); ++n)

{

*drdata = d[n].real();

*didata = d[n].imag();

++drdata;

++didata;

}

return info;

}

};

auto R = Mat<value_type>(N,N);

auto L = Mat<value_type>{};

if(Lr && Li) resize(L,N,N);

auto info = Diag::call(N,M.data(),L.data(),R.data(),dr.data(),di.data());

if(info != 0)

{

//println("M = \n",M);

throw std::runtime_error("Error condition in diagGeneral");

}

struct Unpack

{

void static

call(VectorRef di, MatrixRef Vr, MatrixRef Vi, MatrixRefc V)

{

//Unpack information in V

//back into actual eigenvectors

auto N = di.size();

decltype(N) n = 0;

while(n < N)

{

if(di(n) > 0)

{

//complex eigenvalue pair

column(Vr,n) &= column(V,n);

column(Vr,n+1) &= column(V,n);

column(Vi,n) &= column(V,n+1);

column(Vi,n+1) &= column(V,n+1);

column(Vi,n+1) *= -1;

n += 2;

}

else

{

column(Vr,n) &= column(V,n);

stdx::fill(column(Vi,n),0.);

n += 1;

}

}

}

void static

call(VectorRef di, MatrixRef Vr, MatrixRef Vi, CMatrixRefc V)

{

auto N = di.size();

for(decltype(N) c = 0; c < N; ++c)

for(decltype(N) r = 0; r < N; ++r)

{

Vr(r,c) = V(r,c).real();

Vi(r,c) = V(r,c).imag();

}

}

};

auto Rref = isTransposed(M) ? transpose(R) : makeRef(R);

Unpack::call(makeRef(di),makeRef(Rr),makeRef(Ri),Rref);

if(L)

{

auto Lref = isTransposed(M) ? transpose(L) : makeRef(L);

Unpack::call(makeRef(di),makeRef(Lr),makeRef(Li),Lref);

Error("Inverse step not fully implemented");

//for(auto n : range(N))

// {

// auto facr = column(Lr,n)*column(Rr,n)+column(Li,n)*column(Ri,n);

// auto faci = column(Lr,n)*column(Ri,n)-column(Li,n)*column(Rr,n);

// auto z = Cplx(facr,faci);

// printfln("z %d = %.4E",n,z);

// if(std::abs(z) <= 1E-16) Error("Ill conditioned or non-invertible matrix");

// z = 1./z;

// column(Lr,n) &= column(Lr,n)*z.real()+column(Li,n)*z.imag();

// column(Li,n) &= column(Lr,n)*z.imag()-column(Li,n)*z.real();

// }

}

size_t numpass)

{

auto nkeep = std::min(nrows(M), ncols(M));

auto dots = Vec<V>(nkeep);

for(auto i : range(nkeep))

{

//normalize column i

auto coli = column(M,i);

auto nrm = norm(coli);

if(nrm == 0.0)

{

randomize(coli);

nrm = norm(coli);

}

coli /= nrm;

if(i == 0) continue;

auto Mcols = columns(M,0,i);

auto dotsref = subVector(dots,0,i);

for(auto pass : range1(numpass))

{

// does dotsref &= dag(Mcols) * coli:

auto ccoli = conj(coli);

mult(Mcols,makeRef(ccoli),dotsref,true);

conjugate(dotsref);

// does coli -= Mcols * dotsref:

multSub(Mcols,dotsref,coli);

nrm = norm(coli);

if(nrm < 1E-3) --pass; //orthog is suspect

if(nrm < 1E-10) // What if a subspace was zero in all vectors?

{

randomize(coli);

nrm = norm(coli);

}

coli /= nrm;

}

}

MatrixRefc const& U,

VectorRefc const& D,

MatrixRefc const& V)

{

Matrix Ach(U);

for(auto i : range1(D.size())) column(Ach,i) *= D(i);

Ach = Ach * transpose(V);

Ach -= A;

printfln("relative error with sqrt in low level svd is %.5E",norm(Ach)/norm(A));

Real thresh)

{

auto N = D.size();

if(N <= 1 || thresh <= 0)

{

//println("Got zero thresh");

return make_tuple(true,1);

}

auto D1t = D(0)*thresh;

size_t start = 1;

for(; start < N; ++start)

{

if(D(start) < D1t) break;

}

if(start >= (N-1))

return make_tuple(true,start);

return make_tuple(false,start);

MatRef<T> const& U,

VectorRef const& D,

MatRef<T> const& V,

const Args & args)

{

auto Mr = nrows(M);

auto thresh = args.getReal("SVDThreshold",SVD_THRESH);

//Form 'density matrix' rho

Mat<T> rho, Mconj, tempV, R;

if(isCplx(M))

{

Mconj = conj(M);

rho = M * transpose(Mconj);

}

else

{

rho = M * transpose(M);

}

//Diagonalize rho: evals are squares of singular vals

diagHermitian(rho,U,D);

//Put result of Mt*U==(V*D) in V storage

if(isCplx(M)) mult(transpose(Mconj),U,V);

else mult(transpose(M),U,V);

QR(V, tempV, R, {"Complete=",false, "PositiveDiagonal=",true});

V &= std::move(tempV);

for(decltype(D.size()) i = 0; i < D.size(); ++i)

{

D(i) = std::real(R(i,i));

}

auto [done,start] = checkSVDDone(D,thresh);

if(done) return;

//

//Recursively SVD part of B

//for greater final accuracy

//

auto n = Mr-start;

//reuse rho's storage to avoid allocation

auto mv = move(rho);

reduceCols(mv,n);

auto u = columns(U,start,ncols(U));

auto v = columns(V,start,ncols(V));

//b should be close to diagonal

//but may not be perfect - fix it up below

mult(M,v,mv);

Mat<T> b;

if(isCplx(M)) b = conj(transpose(u))*mv;

else b = transpose(u)*mv;

auto d = subVector(D,start,Mr);

Mat<T> bu(n,n),

bv(n,n);

SVDRefImpl(makeRef(b),makeRef(bu),d,makeRef(bv),args);

//reuse mv's storage to avoid allocation

auto W = move(mv);

mult(u,bu,W);

u &= W;

auto X = v*bv;

v &= X;

MatRef<T> const& U,

VectorRef const& D,

MatRef<T> const& V,

const Args & args)

{

auto Mr = nrows(M), Mc = ncols(M);

if(Mr > Mc)

{

SVDRef(transpose(M),V,D,U,args);

conjugate(V);

conjugate(U);

}

else

{

#ifdef DEBUG

if(!(nrows(U)==Mr && ncols(U)==Mr))

throw std::runtime_error("SVD (ref version), wrong size of U");

if(!(nrows(V)==Mc && ncols(V)==Mr))

throw std::runtime_error("SVD (ref version), wrong size of V");

if(D.size()!=Mr)

throw std::runtime_error("SVD (ref version), wrong size of D");

#endif

auto svdMethod = args.getString("SVDMethod","ITensor");

if(svdMethod=="ITensor")

{

SVDRefImpl(M,U,D,V,args);

}

else if(svdMethod == "gesdd" or svdMethod == "gesvd")

{

SVDRefLAPACK(M,U,D,V,args);

}

else

{

throw std::runtime_error("Unsupported SVD method: "+svdMethod);

}

}

#ifdef CHKSVD

checksvd(M,U,D,V);

#endif

{

LAPACK_INT info = 0;

//

// Scaling: seek ns such that ||F/2^ns|| < 1/2

// and set scale = 1/2^ns

//

auto ns = dlange_wrapper('I',N,N,F.data());// infinite norm of the matrix to be exponentiated

#ifdef DEBUG

if(ns == 0) throw std::runtime_error("padeExp: null input matrix");

#endif

ns = std::max(0,(int)(std::log(ns)/log(2))+2);

Real scale = std::pow(2,-ns);

Real scale2 = scale*scale;

//

// Compute Pade coefficient

//

std::vector<Real> coef(ideg+1);

coef[0] = 1.0;

for(int k = 1; k <= ideg; ++k)

{

coef[k] = coef[k-1]*((double)(ideg+1-k)/(double)(k*(2*ideg+1-k)));

}

//

// H^2 = scale2*F*F

//

auto F2 = Mat<Real>(N,N);

gemm(F,F,makeRef(F2),scale2,0);

//

// Initialize P and Q

//

auto P = Mat<Real>(N,N);

auto Q = Mat<Real>(N,N);

for(auto j : range(N))

{

Q(j,j) = coef[ideg];

P(j,j) = coef[ideg-1];

}

//

// Horner evaluation of the irreducible fraction:

// Apply Horner rule

//

bool odd = true;

for(int k = ideg-1; k > 0; --k)

{

if(odd)

{

Q = Q*F2;

for(auto j : range(N))

{

Q(j,j) += coef[k-1];

}

}

else

{

P = P*F2;

for(auto j : range(N))

{

P(j,j) += coef[k-1];

}

}

odd = !odd;

}

//

// Horner evaluation of the irreducible fraction:

// Obtain (+/-)(I+2*(P\Q))

//

if(odd)

{

Q = scale*Q*F;

}

else

{

P = scale*P*F;

}

Q -= P;

info = dgesv_wrapper(N,N,Q.data(),P.data());

if(info != 0) return info;

P *= 2.0;

for(auto j : range(N))

{

P(j,j) += 1.0;

}

if(ns == 0 && odd)

{

P *= -1.0;

}

//

// Squaring: exp(F) = (exp(F))^(2^ns)

//

for(int k = 1; k <= ns; ++k)

{

P = P*P;

}

//auto pend = P.data()+P.size();

//auto f = Fdata;

//for(auto p = P.data(); p != pend; ++p,++f)

// {

// *f = *p;

// }

F &= P;//deep copy

return info;

{

LAPACK_INT info = 0;

//

// Scaling: seek ns such that ||F/2^ns|| < 1/2

// and set scale = 1/2^ns

//

auto ns = zlange_wrapper('I',N,N,F.data());

#ifdef DEBUG

if(ns == 0) throw std::runtime_error("padeExp: null input matrix");

#endif

ns = std::max(0,(int)(std::log(ns)/log(2))+2);

Real scale = std::pow(2,-ns);

Real scale2 = scale*scale;

//

// Compute Pade coefficient

//

std::vector<Real> coef(ideg+1);

coef[0] = 1.0;

for(int k = 1; k <= ideg; ++k)

{

coef[k] = coef[k-1]*((double)(ideg+1-k)/(double)(k*(2*ideg+1-k)));

}

//

// H^2 = scale2*F*F

//

auto F2 = Mat<Cplx>(N,N);

gemm(F,F,makeRef(F2),scale2,0);

//

// Initialize P and Q

//

auto P = Mat<Cplx>(N,N);

auto Q = Mat<Cplx>(N,N);

for(auto j : range(N))

{

Q(j,j) = coef[ideg];

P(j,j) = coef[ideg-1];

}

//

// Horner evaluation of the irreducible fraction:

// Apply Horner rule

//

bool odd = true;

for(int k = ideg-1; k > 0; --k)

{

if(odd)

{

Q = Q*F2;

for(auto j : range(N))

{

Q(j,j) += coef[k-1];

}

}

else

{

P = P*F2;

for(auto j : range(N))

{

P(j,j) += coef[k-1];

}

}

odd = !odd;

}