本文整理汇总了C++中DiagonalMatrix::element方法的典型用法代码示例。如果您正苦于以下问题：C++ DiagonalMatrix::element方法的具体用法？C++ DiagonalMatrix::element怎么用？C++ DiagonalMatrix::element使用的例子？那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类DiagonalMatrix的用法示例。
在下文中一共展示了DiagonalMatrix::element方法的5个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: _sortEigens
void PrincipalComponentsAnalysis::_sortEigens(Matrix& eigenVectors,
DiagonalMatrix& eigenValues)
{
// simple bubble sort, slow, but I don't expect very large matrices. Lots of room for
// improvement.
bool change = true;
while (change)
{
change = false;
for (int c = 0; c < eigenVectors.Ncols() - 1; c++)
{
if (eigenValues.element(c) < eigenValues.element(c + 1))
{
std::swap(eigenValues.element(c), eigenValues.element(c + 1));
for (int r = 0; r < eigenVectors.Nrows(); r++)
{
std::swap(eigenVectors.element(r, c), eigenVectors.element(r, c + 1));
}
}
}
}
}示例2: Jacobi
void Jacobi(const SymmetricMatrix& X, DiagonalMatrix& D, SymmetricMatrix& A,
Matrix& V, bool eivec)
{
Real epsilon = FloatingPointPrecision::Epsilon();
Tracer et("Jacobi");
REPORT
int n = X.Nrows(); DiagonalMatrix B(n), Z(n); D.resize(n); A = X;
if (eivec) { REPORT V.resize(n,n); D = 1.0; V = D; }
B << A; D = B; Z = 0.0; A.Inject(Z);
bool converged = false;
for (int i=1; i<=50; i++)
{
Real sm=0.0; Real* a = A.Store(); int p = A.Storage();
while (p--) sm += fabs(*a++); // have previously zeroed diags
if (sm==0.0) { REPORT converged = true; break; }
Real tresh = (i<4) ? 0.2 * sm / square(n) : 0.0; a = A.Store();
for (p = 0; p < n; p++)
{
Real* ap1 = a + (p*(p+1))/2;
Real& zp = Z.element(p); Real& dp = D.element(p);
for (int q = p+1; q < n; q++)
{
Real* ap = ap1; Real* aq = a + (q*(q+1))/2;
Real& zq = Z.element(q); Real& dq = D.element(q);
Real& apq = A.element(q,p);
Real g = 100 * fabs(apq); Real adp = fabs(dp); Real adq = fabs(dq);
if (i>4 && g < epsilon*adp && g < epsilon*adq) { REPORT apq = 0.0; }
else if (fabs(apq) > tresh)
{
REPORT
Real t; Real h = dq - dp; Real ah = fabs(h);
if (g < epsilon*ah) { REPORT t = apq / h; }
else
{
REPORT
Real theta = 0.5 * h / apq;
t = 1.0 / ( fabs(theta) + sqrt(1.0 + square(theta)) );
if (theta<0.0) { REPORT t = -t; }
}
Real c = 1.0 / sqrt(1.0 + square(t)); Real s = t * c;
Real tau = s / (1.0 + c); h = t * apq;
zp -= h; zq += h; dp -= h; dq += h; apq = 0.0;
int j = p;
while (j--)
{
g = *ap; h = *aq;
*ap++ = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
}
int ip = p+1; j = q-ip; ap += ip++; aq++;
while (j--)
{
g = *ap; h = *aq;
*ap = g-s*(h+g*tau); *aq++ = h+s*(g-h*tau);
ap += ip++;
}
if (q < n-1) // last loop is non-empty
{
int iq = q+1; j = n-iq; ap += ip++; aq += iq++;
for (;;)
{
g = *ap; h = *aq;
*ap = g-s*(h+g*tau); *aq = h+s*(g-h*tau);
if (!(--j)) break;
ap += ip++; aq += iq++;
}
}
if (eivec)
{
REPORT
RectMatrixCol VP(V,p); RectMatrixCol VQ(V,q);
Rotate(VP, VQ, tau, s);
}
}
}
}
B = B + Z; D = B; Z = 0.0;
}
if (!converged) Throw(ConvergenceException(X));
if (eivec) SortSV(D, V, true);
else SortAscending(D);
}示例3: subspace_h
void SpinAdapted::Linear::block_davidson(vector<Wavefunction>& b, DiagonalMatrix& h_diag, double normtol, const bool &warmUp, Davidson_functor& h_multiply, bool& useprecond, int currentRoot, std::vector<Wavefunction> &lowerStates)
{
pout.precision (12);
#ifndef SERIAL
mpi::communicator world;
#endif
int iter = 0;
double levelshift = 0.0;
int nroots = b.size();
//normalise all the guess roots
if(mpigetrank() == 0) {
for(int i=0; i<nroots; ++i)
{
for(int j=0; j<i; ++j)
{
double overlap = DotProduct(b[j], b[i]);
ScaleAdd(-overlap, b[j], b[i]);
}
Normalise(b[i]);
}
//if we are doing state specific, lowerstates has lower energy states
if (lowerStates.size() != 0) {
for (int i=0; i<lowerStates.size(); i++) {
double overlap = DotProduct(b[0], lowerStates[i]);
ScaleAdd(-overlap/DotProduct(lowerStates[i], lowerStates[i]), lowerStates[i], b[0]);
}
Normalise(b[0]);
}
}
vector<Wavefunction> sigma;
int converged_roots = 0;
int maxiter = h_diag.Ncols() - lowerStates.size();
while(true)
{
if (dmrginp.outputlevel() > 0)
pout << "\t\t\t Davidson Iteration :: " << iter << endl;
++iter;
dmrginp.hmultiply -> start();
int sigmasize, bsize;
if (mpigetrank() == 0) {
sigmasize = sigma.size();
bsize = b.size();
}
#ifndef SERIAL
mpi::broadcast(world, sigmasize, 0);
mpi::broadcast(world, bsize, 0);
#endif
//multiply all guess vectors with hamiltonian c = Hv
for(int i=sigmasize; i<bsize; ++i) {
Wavefunction sigmai, bi;
Wavefunction* sigmaptr=&sigmai, *bptr = &bi;
if (mpigetrank() == 0) {
sigma.push_back(b[i]);
sigma[i].Clear();
sigmaptr = &sigma[i];
bptr = &b[i];
}
#ifndef SERIAL
mpi::broadcast(world, *bptr, 0);
#endif
if (mpigetrank() != 0) {
sigmai = bi;
sigmai.Clear();
}
h_multiply(*bptr, *sigmaptr);
//if (mpigetrank() == 0) {
// cout << *bptr << endl;
// cout << *sigmaptr << endl;
//}
}
dmrginp.hmultiply -> stop();
Wavefunction r;
DiagonalMatrix subspace_eigenvalues;
if (mpigetrank() == 0) {
Matrix subspace_h(b.size(), b.size());
for (int i = 0; i < b.size(); ++i)
for (int j = 0; j <= i; ++j) {
subspace_h.element(i, j) = DotProduct(b[i], sigma[j]);
subspace_h.element(j, i) = subspace_h.element(i, j);
}
Matrix alpha;
diagonalise(subspace_h, subspace_eigenvalues, alpha);
if (dmrginp.outputlevel() > 0) {
for (int i = 1; i <= subspace_eigenvalues.Ncols (); ++i)
pout << "\t\t\t " << i << " :: " << subspace_eigenvalues(i,i)+dmrginp.get_coreenergy() << endl;
//.........这里部分代码省略.........
示例4: make_implicit_svd
//==========================================================================
void make_implicit_svd(vector<vector<double> >& mat,
vector<double>& b, double& sigma_min)
//==========================================================================
{
int rows = (int)mat.size();
int cols = (int)mat[0].size();
// cout << "Rows = " << rows << endl
// << "Cols = " << cols << endl;
Matrix nmat;
nmat.ReSize(rows, cols);
for (int i = 0; i < rows; ++i) {
for (int j = 0; j < cols; ++j) {
nmat.element(i, j) = mat[i][j];
}
}
// Check if mat has enough rows. If not fill out with zeros.
if (rows < cols) {
RowVector zero(cols);
zero = 0.0; // Initializes zero to a null-vector.
for (int i = rows; i < cols; ++i) {
nmat &= zero; // & means horizontal concatenation in newmat
}
}
// Perform SVD.
// cout << "Running SVD..." << endl;
static DiagonalMatrix diag;
static Matrix V;
Try {
SVD(nmat, diag, nmat, V);
} CatchAll {
cout << Exception::what() << endl;
b = vector<double>(cols, 0.0);
sigma_min = -1.0;
return;
}
// // Write out singular values.
// cout << "Singular values:" << endl;
// for (int ik = 0; ik < cols; ik++)
// cout << ik << "\t" << diag.element(ik, ik) << endl;
// // Write out info about singular values
// double s_min = diag.element(cols-1, cols-1);
// double s_max = diag.element(0, 0);
// cout << "Implicitization:" << endl
// << "s_min = " << s_min << endl
// << "s_max = " << s_max << endl
// << "Ratio of s_min/s_max = " << s_min/s_max << endl;
// // Find square sum of singular values
// double sum = 0.0;
// for (int i = 0; i < cols; i++)
// sum += diag.element(i, i) * diag.element(i, i);
// sum = sqrt(sum);
// cout << "Square sum = " << sum << endl;
// Get the appropriate null-vector and corresponding singular value
const double eps = 1.0e-15;
double tol = cols * fabs(diag.element(0, 0)) * eps;
int nullvec = 0;
for (int i = 0; i < cols-1; ++i) {
if (fabs(diag.element(i, i)) > tol) {
++nullvec;
}
}
sigma_min = diag.element(nullvec, nullvec);
// cout << "Null-vector: " << nullvec << endl
// << "sigma_min = " << sigma_min << endl;
// Set the coefficients
b.resize(cols);
for (int jk = 0; jk < cols; ++jk)
b[jk] = V.element(jk, nullvec);
return;
}示例5: et
static void tql2(DiagonalMatrix& D, DiagonalMatrix& E, Matrix& Z)
{
Tracer et("Evalue(tql2)");
Real eps = FloatingPointPrecision::Epsilon();
int n = D.Nrows(); Real* z = Z.Store(); int l;
for (l=1; l<n; l++) E.element(l-1) = E.element(l);
Real b = 0.0; Real f = 0.0; E.element(n-1) = 0.0;
for (l=0; l<n; l++)
{
int i,j;
Real& dl = D.element(l); Real& el = E.element(l);
Real h = eps * ( fabs(dl) + fabs(el) );
if (b < h) b = h;
int m;
for (m=l; m<n; m++) if (fabs(E.element(m)) <= b) break;
bool test = false;
for (j=0; j<30; j++)
{
if (m==l) { test = true; break; }
Real& dl1 = D.element(l+1);
Real g = dl; Real p = (dl1-g) / (2.0*el); Real r = sqrt(p*p + 1.0);
dl = el / (p < 0.0 ? p-r : p+r); Real h = g - dl; f += h;
Real* dlx = &dl1; i = n-l-1; while (i--) *dlx++ -= h;
p = D.element(m); Real c = 1.0; Real s = 0.0;
for (i=m-1; i>=l; i--)
{
Real ei = E.element(i); Real di = D.element(i);
Real& ei1 = E.element(i+1);
g = c * ei; h = c * p;
if ( fabs(p) >= fabs(ei))
{
c = ei / p; r = sqrt(c*c + 1.0);
ei1 = s*p*r; s = c/r; c = 1.0/r;
}
else
{
c = p / ei; r = sqrt(c*c + 1.0);
ei1 = s * ei * r; s = 1.0/r; c /= r;
}
p = c * di - s*g; D.element(i+1) = h + s * (c*g + s*di);
Real* zki = z + i; Real* zki1 = zki + 1; int k = n;
while (k--)
{
h = *zki1; *zki1 = s*(*zki) + c*h; *zki = c*(*zki) - s*h;
zki += n; zki1 += n;
}
}
el = s*p; dl = c*p;
if (fabs(el) <= b) { test = true; break; }
}
if (!test) Throw ( ConvergenceException(D) );
dl += f;
}
for (int i=0; i<n; i++)
{
int k = i; Real p = D.element(i);
for (int j=i+1; j<n; j++)
{ if (D.element(j) < p) { k = j; p = D.element(j); } }
if (k != i)
{
D.element(k) = D.element(i); D.element(i) = p; int j = n;
Real* zji = z + i; Real* zjk = z + k;
while (j--) { p = *zji; *zji = *zjk; *zjk = p; zji += n; zjk += n; }
}
}
}