* @note Requires LAPACK support.
*/
bool compute_all_condition_numbers;
-
+
/**
* Compute all eigenvalues of
* the projected matrix.
* @note Requires LAPACK support.
*/
bool compute_eigenvalues;
-
+
/**
* Constructor. Initialize data
* fields. Confer the description of
VECTOR *Vp;
VECTOR *Vz;
VECTOR *VAp;
-
+
/**
* Within the iteration loop, the
* square of the residual vector is
compute_all_condition_numbers(compute_all_condition_numbers),
compute_eigenvalues(compute_eigenvalues)
{}
-
+
template <class VECTOR>
SolverControl::State conv=SolverControl::iterate;
deallog.push("cg");
-
+
// Memory allocation
Vr = this->memory.alloc();
Vp = this->memory.alloc();
| additional_data.compute_all_condition_numbers
| additional_data.compute_eigenvalues;
double eigen_beta_alpha = 0;
-
+
// vectors used for eigenvalue
// computations
std::vector<double> diagonal;
std::vector<double> offdiagonal;
-
+
try {
// define some aliases for simpler access
- VECTOR& g = *Vr;
- VECTOR& h = *Vp;
- VECTOR& d = *Vz;
+ VECTOR& g = *Vr;
+ VECTOR& h = *Vp;
+ VECTOR& d = *Vz;
VECTOR& Ad = *VAp;
// resize the vectors, but do not set
// the values since they'd be overwritten
// library
int it=0;
double res,gh,alpha,beta;
-
+
// compute residual. if vector is
// zero, then short-circuit the
// full computation
res = g.l2_norm();
conv = this->control().check(0,res);
- if (conv)
+ if (conv)
{
cleanup();
return;
}
precondition.vmult(h,g);
-
+
d.equ(-1.,h);
-
+
gh = g*h;
-
+
while (conv == SolverControl::iterate)
{
it++;
A.vmult(Ad,d);
-
+
alpha = d*Ad;
alpha = gh/alpha;
-
+
g.add(alpha,Ad);
x.add(alpha,d );
res = g.l2_norm();
-
+
print_vectors(it, x, g, d);
-
+
conv = this->control().check(it,res);
if (conv)
break;
precondition.vmult(h,g);
-
+
beta = gh;
gh = g*h;
beta = gh/beta;
-
+
if (additional_data.log_coefficients)
deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
// set up the vectors
deallog << "Condition number estimate: " <<
T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
}
-
+
d.sadd(beta,-1.,h);
}
}
* Set the element (<i>i,j</i>)
* to @p value.
*
- * This function is able to insert
- * new elements into the matrix as
- * long as compress() has not been
- * called, so the sparsity pattern
- * will be extended. When compress()
- * is called for the first time, then
- * this is no longer possible and an
- * insertion of elements at positions
- * which have not been initialized
- * will throw an exception. Moreover,
- * if <tt>value</tt> is not a finite
- * number an exception is thrown.
+ * This function is able to insert new
+ * elements into the matrix as long as
+ * compress() has not been called, so
+ * the sparsity pattern will be
+ * extended. When compress() is called
+ * for the first time, then this is no
+ * longer possible and an insertion of
+ * elements at positions which have not
+ * been initialized will throw an
+ * exception. Note that in case
+ * elements need to be inserted, it is
+ * mandatory that elements are inserted
+ * only once. Otherwise, the elements
+ * will actually be added in the end
+ * (since it is not possible to
+ * efficiently find values to the same
+ * entry before compress() has been
+ * called). In the case that an element
+ * is set more than once, initialize
+ * the matrix with a sparsity pattern
+ * first.
*/
void set (const unsigned int i,
const unsigned int j,