ROL
ROL_Sketch.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_SKETCH_H
45#define ROL_SKETCH_H
46
47#include "ROL_Vector.hpp"
48#include "ROL_LinearAlgebra.hpp"
49#include "ROL_LAPACK.hpp"
50#include "ROL_Types.hpp"
51#include <random>
52#include <chrono>
53
61namespace ROL {
62
63template <class Real>
64class Sketch {
65private:
66 std::vector<Ptr<Vector<Real>>> Upsilon_, Phi_, Y_;
67 LA::Matrix<Real> Omega_, Psi_, X_, Z_, C_;
68
70
71 const Real orthTol_;
72 const int orthIt_;
73
74 const bool truncate_;
75
76 LAPACK<int,Real> lapack_;
77
79
80 Ptr<std::ostream> out_;
81
82 Ptr<Elementwise::NormalRandom<Real>> nrand_;
83 Ptr<std::mt19937_64> gen_;
84 Ptr<std::normal_distribution<Real>> dist_;
85
86 int computeP(void) {
87 int INFO(0);
88 if (!flagP_) {
89 // Solve least squares problem using LAPACK
90 int M = ncol_;
91 int N = k_;
92 int K = std::min(M,N);
93 int LDA = M;
94 std::vector<Real> TAU(K);
95 std::vector<Real> WORK(1);
96 int LWORK = -1;
97 // Compute QR factorization of X
98 lapack_.GEQRF(M,N,X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
99 LWORK = static_cast<int>(WORK[0]);
100 WORK.resize(LWORK);
101 lapack_.GEQRF(M,N,X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
102 // Generate Q
103 LWORK = -1;
104 lapack_.ORGQR(M,N,K,X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
105 LWORK = static_cast<int>(WORK[0]);
106 WORK.resize(LWORK);
107 lapack_.ORGQR(M,N,K,X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
108 flagP_ = true;
109 }
110 return INFO;
111 }
112
113 void mgs2(std::vector<Ptr<Vector<Real>>> &Y) const {
114 const Real one(1);
115 Real rjj(0), rij(0);
116 std::vector<Real> normQ(k_,0);
117 bool flag(true);
118 for (int j = 0; j < k_; ++j) {
119 rjj = Y[j]->norm();
120 if (rjj > ROL_EPSILON<Real>()) { // Ignore update if Y[i] is zero.
121 for (int k = 0; k < orthIt_; ++k) {
122 for (int i = 0; i < j; ++i) {
123 rij = Y[i]->dot(*Y[j]);
124 Y[j]->axpy(-rij,*Y[i]);
125 }
126 normQ[j] = Y[j]->norm();
127 flag = true;
128 for (int i = 0; i < j; ++i) {
129 rij = std::abs(Y[i]->dot(*Y[j]));
130 if (rij > orthTol_*normQ[j]*normQ[i]) {
131 flag = false;
132 break;
133 }
134 }
135 if (flag) {
136 break;
137 }
138 }
139 }
140 rjj = normQ[j];
141 Y[j]->scale(one/rjj);
142 }
143 }
144
145 int computeQ(void) {
146 if (!flagQ_) {
147 mgs2(Y_);
148 flagQ_ = true;
149 }
150 return 0;
151 }
152
153 int LSsolver(LA::Matrix<Real> &A, LA::Matrix<Real> &B, const bool trans = false) const {
154 int flag(0);
155 char TRANS = (trans ? 'T' : 'N');
156 int M = A.numRows();
157 int N = A.numCols();
158 int NRHS = B.numCols();
159 int LDA = M;
160 int LDB = std::max(M,N);
161 std::vector<Real> WORK(1);
162 int LWORK = -1;
163 int INFO;
164 lapack_.GELS(TRANS,M,N,NRHS,A.values(),LDA,B.values(),LDB,&WORK[0],LWORK,&INFO);
165 flag += INFO;
166 LWORK = static_cast<int>(WORK[0]);
167 WORK.resize(LWORK);
168 lapack_.GELS(TRANS,M,N,NRHS,A.values(),LDA,B.values(),LDB,&WORK[0],LWORK,&INFO);
169 flag += INFO;
170 return flag;
171 }
172
173 int lowRankApprox(LA::Matrix<Real> &A, const int r) const {
174 const Real zero(0);
175 char JOBU = 'S';
176 char JOBVT = 'S';
177 int M = A.numRows();
178 int N = A.numCols();
179 int K = std::min(M,N);
180 int LDA = M;
181 std::vector<Real> S(K);
182 LA::Matrix<Real> U(M,K);
183 int LDU = M;
184 LA::Matrix<Real> VT(K,N);
185 int LDVT = K;
186 std::vector<Real> WORK(1), WORK0(1);
187 int LWORK = -1;
188 int INFO;
189 lapack_.GESVD(JOBU,JOBVT,M,N,A.values(),LDA,&S[0],U.values(),LDU,VT.values(),LDVT,&WORK[0],LWORK,&WORK0[0],&INFO);
190 LWORK = static_cast<int>(WORK[0]);
191 WORK.resize(LWORK);
192 lapack_.GESVD(JOBU,JOBVT,M,N,A.values(),LDA,&S[0],U.values(),LDU,VT.values(),LDVT,&WORK[0],LWORK,&WORK0[0],&INFO);
193 for (int i = 0; i < M; ++i) {
194 for (int j = 0; j < N; ++j) {
195 A(i,j) = zero;
196 for (int k = 0; k < r; ++k) {
197 A(i,j) += S[k] * U(i,k) * VT(k,j);
198 }
199 }
200 }
201 return INFO;
202 }
203
204 int computeC(void) {
205 int infoP(0), infoQ(0), infoLS1(0), infoLS2(0), infoLRA(0);
206 infoP = computeP();
207 infoQ = computeQ();
208 if (!flagC_) {
209 const Real zero(0);
210 LA::Matrix<Real> L(s_,k_), R(s_,k_);
211 for (int i = 0; i < s_; ++i) {
212 for (int j = 0; j < k_; ++j) {
213 L(i,j) = Phi_[i]->dot(*Y_[j]);
214 R(i,j) = zero;
215 for (int k = 0; k < ncol_; ++k) {
216 R(i,j) += Psi_(k,i) * X_(k,j);
217 }
218 }
219 }
220 // Solve least squares problems using LAPACK
221 infoLS1 = LSsolver(L,Z_,false);
222 LA::Matrix<Real> Z(s_,k_);
223 for (int i = 0; i < k_; ++i) {
224 for (int j = 0; j < s_; ++j) {
225 Z(j,i) = Z_(i,j);
226 }
227 }
228 infoLS2 = LSsolver(R,Z,false);
229 for (int i = 0; i < k_; ++i) {
230 for (int j = 0; j < k_; ++j) {
231 C_(j,i) = Z(i,j);
232 }
233 }
234 // Compute best rank r approximation
235 if (truncate_) {
236 infoLRA = lowRankApprox(C_,rank_);
237 }
238 // Set flag
239 flagC_ = true;
240 }
241 return std::abs(infoP)+std::abs(infoQ)+std::abs(infoLS1)
242 +std::abs(infoLS2)+std::abs(infoLRA);
243 }
244
245 void reset(void) {
246 const Real zero(0);
247 // Randomize Upsilon, Omega, Psi and Phi, and zero X, Y and Z
248 X_.scale(zero); Z_.scale(zero); C_.scale(zero);
249 for (int i = 0; i < s_; ++i) {
250 Phi_[i]->applyUnary(*nrand_);
251 for (int j = 0; j < ncol_; ++j) {
252 Psi_(j,i) = (*dist_)(*gen_);
253 }
254 }
255 for (int i = 0; i < k_; ++i) {
256 Y_[i]->zero();
257 Upsilon_[i]->applyUnary(*nrand_);
258 for (int j = 0; j < ncol_; ++j) {
259 Omega_(j,i) = (*dist_)(*gen_);
260 }
261 }
262 }
263
264// void reset(void) {
265// const Real zero(0), one(1), a(2), b(-1);
266// Real x(0);
267// // Randomize Upsilon, Omega, Psi and Phi, and zero X, Y and Z
268// X_.scale(zero); Z_.scale(zero); C_.scale(zero);
269// for (int i = 0; i < s_; ++i) {
270// Phi_[i]->randomize(-one,one);
271// for (int j = 0; j < ncol_; ++j) {
272// x = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
273// Psi_(j,i) = a*x + b;
274// }
275// }
276// for (int i = 0; i < k_; ++i) {
277// Y_[i]->zero();
278// Upsilon_[i]->randomize(-one,one);
279// for (int j = 0; j < ncol_; ++j) {
280// x = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
281// Omega_(j,i) = a*x + b;
282// }
283// }
284// }
285
286public:
287 virtual ~Sketch(void) {}
288
289 Sketch(const Vector<Real> &x, const int ncol, const int rank,
290 const Real orthTol = 1e-8, const int orthIt = 2,
291 const bool truncate = false,
292 const unsigned dom_seed = 0, const unsigned rng_seed = 0)
293 : ncol_(ncol), orthTol_(orthTol), orthIt_(orthIt),
294 truncate_(truncate), flagP_(false), flagQ_(false), flagC_(false),
295 out_(nullPtr) {
296 Real mu(0), sig(1);
297 nrand_ = makePtr<Elementwise::NormalRandom<Real>>(mu,sig,dom_seed);
298 unsigned seed = rng_seed;
299 if (seed == 0) {
300 seed = std::chrono::system_clock::now().time_since_epoch().count();
301 }
302 gen_ = makePtr<std::mt19937_64>(seed);
303 dist_ = makePtr<std::normal_distribution<Real>>(mu,sig);
304 // Compute reduced dimensions
305 maxRank_ = std::min(ncol_, x.dimension());
306 rank_ = std::min(rank, maxRank_);
307 k_ = std::min(2*rank_+1, maxRank_);
308 s_ = std::min(2*k_ +1, maxRank_);
309 // Initialize matrix storage
310 Upsilon_.clear(); Phi_.clear(); Y_.clear();
311 Omega_.reshape(ncol_,k_); Psi_.reshape(ncol_,s_);
312 X_.reshape(ncol_,k_); Z_.reshape(s_,s_); C_.reshape(k_,k_);
313 for (int i = 0; i < k_; ++i) {
314 Y_.push_back(x.clone());
315 Upsilon_.push_back(x.clone());
316 }
317 for (int i = 0; i < s_; ++i) {
318 Phi_.push_back(x.clone());
319 }
320 // Randomize Psi and Omega and zero W and Y
321 reset();
322 }
323
324 void setStream(Ptr<std::ostream> &out) {
325 out_ = out;
326 }
327
328 void setRank(const int rank) {
329 rank_ = std::min(rank, maxRank_);
330 // Compute reduced dimensions
331 int sold = s_, kold = k_;
332 k_ = std::min(2*rank_+1, maxRank_);
333 s_ = std::min(2*k_ +1, maxRank_);
334 Omega_.reshape(ncol_,k_); Psi_.reshape(ncol_,s_);
335 X_.reshape(ncol_,k_); Z_.reshape(s_,s_); C_.reshape(k_,k_);
336 if (s_ > sold) {
337 for (int i = sold; i < s_; ++i) {
338 Phi_.push_back(Phi_[0]->clone());
339 }
340 }
341 if (k_ > kold) {
342 for (int i = kold; i < k_; ++i) {
343 Y_.push_back(Y_[0]->clone());
344 Upsilon_.push_back(Upsilon_[0]->clone());
345 }
346 }
347 update();
348 if ( out_ != nullPtr ) {
349 *out_ << std::string(80,'=') << std::endl;
350 *out_ << " ROL::Sketch::setRank" << std::endl;
351 *out_ << " **** Rank = " << rank_ << std::endl;
352 *out_ << " **** k = " << k_ << std::endl;
353 *out_ << " **** s = " << s_ << std::endl;
354 *out_ << std::string(80,'=') << std::endl;
355 }
356 }
357
358 void update(void) {
359 flagP_ = false;
360 flagQ_ = false;
361 flagC_ = false;
362 // Randomize Psi and Omega and zero W and Y
363 reset();
364 }
365
366 int advance(const Real nu, Vector<Real> &h, const int col, const Real eta = 1.0) {
367 // Check to see if col is less than ncol_
368 if ( col >= ncol_ || col < 0 ) {
369 // Input column index out of range!
370 return 1;
371 }
372 if (!flagP_ && !flagQ_ && !flagC_) {
373 for (int i = 0; i < k_; ++i) {
374 // Update X
375 for (int j = 0; j < ncol_; ++j) {
376 X_(j,i) *= eta;
377 }
378 X_(col,i) += nu*h.dot(*Upsilon_[i]);
379 // Update Y
380 Y_[i]->scale(eta);
381 Y_[i]->axpy(nu*Omega_(col,i),h);
382 }
383 // Update Z
384 Real hphi(0);
385 for (int i = 0; i < s_; ++i) {
386 hphi = h.dot(*Phi_[i]);
387 for (int j = 0; j < s_; ++j) {
388 Z_(i,j) *= eta;
389 Z_(i,j) += nu*Psi_(col,j)*hphi;
390 }
391 }
392 if ( out_ != nullPtr ) {
393 *out_ << std::string(80,'=') << std::endl;
394 *out_ << " ROL::Sketch::advance" << std::endl;
395 *out_ << " **** col = " << col << std::endl;
396 *out_ << " **** norm(h) = " << h.norm() << std::endl;
397 *out_ << std::string(80,'=') << std::endl;
398 }
399 }
400 else {
401 // Reconstruct has already been called!
402 return 1;
403 }
404 return 0;
405 }
406
407 int reconstruct(Vector<Real> &a, const int col) {
408 // Check to see if col is less than ncol_
409 if ( col >= ncol_ || col < 0 ) {
410 // Input column index out of range!
411 return 2;
412 }
413 const Real zero(0);
414 int flag(0);
415 // Compute QR factorization of X store in X
416 flag = computeP();
417 if (flag > 0 ) {
418 return 3;
419 }
420 // Compute QR factorization of Y store in Y
421 flag = computeQ();
422 if (flag > 0 ) {
423 return 4;
424 }
425 // Compute (Phi Q)\Z/(Psi P)* store in C
426 flag = computeC();
427 if (flag > 0 ) {
428 return 5;
429 }
430 // Recover sketch
431 a.zero();
432 Real coeff(0);
433 for (int i = 0; i < k_; ++i) {
434 coeff = zero;
435 for (int j = 0; j < k_; ++j) {
436 coeff += C_(i,j) * X_(col,j);
437 }
438 a.axpy(coeff,*Y_[i]);
439 }
440 if ( out_ != nullPtr ) {
441 *out_ << std::string(80,'=') << std::endl;
442 *out_ << " ROL::Sketch::reconstruct" << std::endl;
443 *out_ << " **** col = " << col << std::endl;
444 *out_ << " **** norm(a) = " << a.norm() << std::endl;
445 *out_ << std::string(80,'=') << std::endl;
446 }
447 return 0;
448 }
449
450 bool test(const int rank, std::ostream &outStream = std::cout, const int verbosity = 0) {
451 const Real one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
452 // Initialize low rank factors
453 std::vector<Ptr<Vector<Real>>> U(rank);
454 LA::Matrix<Real> V(ncol_,rank);
455 for (int i = 0; i < rank; ++i) {
456 U[i] = Y_[0]->clone();
457 U[i]->randomize(-one,one);
458 for (int j = 0; j < ncol_; ++j) {
459 V(j,i) = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
460 }
461 }
462 // Initialize A and build sketch
463 update();
464 std::vector<Ptr<Vector<Real>>> A(ncol_);
465 for (int i = 0; i < ncol_; ++i) {
466 A[i] = Y_[0]->clone(); A[i]->zero();
467 for (int j = 0; j < rank; ++j) {
468 A[i]->axpy(V(i,j),*U[j]);
469 }
470 advance(one,*A[i],i,one);
471 }
472 // Test QR decomposition of X
473 bool flagP = testP(outStream, verbosity);
474 // Test QR decomposition of Y
475 bool flagQ = testQ(outStream, verbosity);
476 // Test reconstruction of A
477 Real nerr(0), maxerr(0);
478 Ptr<Vector<Real>> err = Y_[0]->clone();
479 for (int i = 0; i < ncol_; ++i) {
480 reconstruct(*err,i);
481 err->axpy(-one,*A[i]);
482 nerr = err->norm();
483 maxerr = (nerr > maxerr ? nerr : maxerr);
484 }
485 if (verbosity > 0) {
486 std::ios_base::fmtflags oflags(outStream.flags());
487 outStream << std::scientific << std::setprecision(3) << std::endl;
488 outStream << " TEST RECONSTRUCTION: Max Error = "
489 << std::setw(12) << std::right << maxerr
490 << std::endl << std::endl;
491 outStream.flags(oflags);
492 }
493 return flagP & flagQ & (maxerr < tol ? true : false);
494 }
495
496private:
497
498 // Test functions
499 bool testQ(std::ostream &outStream = std::cout, const int verbosity = 0) {
500 const Real one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
501 computeQ();
502 Real qij(0), err(0), maxerr(0);
503 std::ios_base::fmtflags oflags(outStream.flags());
504 if (verbosity > 0) {
505 outStream << std::scientific << std::setprecision(3);
506 }
507 if (verbosity > 1) {
508 outStream << std::endl
509 << " Printing Q'Q...This should be approximately equal to I"
510 << std::endl << std::endl;
511 }
512 for (int i = 0; i < k_; ++i) {
513 for (int j = 0; j < k_; ++j) {
514 qij = Y_[i]->dot(*Y_[j]);
515 err = (i==j ? std::abs(qij-one) : std::abs(qij));
516 maxerr = (err > maxerr ? err : maxerr);
517 if (verbosity > 1) {
518 outStream << std::setw(12) << std::right << qij;
519 }
520 }
521 if (verbosity > 1) {
522 outStream << std::endl;
523 }
524 if (maxerr > tol) {
525 break;
526 }
527 }
528 if (verbosity > 0) {
529 outStream << std::endl << " TEST ORTHOGONALIZATION: Max Error = "
530 << std::setw(12) << std::right << maxerr
531 << std::endl;
532 outStream.flags(oflags);
533 }
534 return (maxerr < tol ? true : false);
535 }
536
537 bool testP(std::ostream &outStream = std::cout, const int verbosity = 0) {
538 const Real zero(0), one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
539 computeP();
540 Real qij(0), err(0), maxerr(0);
541 std::ios_base::fmtflags oflags(outStream.flags());
542 if (verbosity > 0) {
543 outStream << std::scientific << std::setprecision(3);
544 }
545 if (verbosity > 1) {
546 outStream << std::endl
547 << " Printing P'P...This should be approximately equal to I"
548 << std::endl << std::endl;
549 }
550 for (int i = 0; i < k_; ++i) {
551 for (int j = 0; j < k_; ++j) {
552 qij = zero;
553 for (int k = 0; k < ncol_; ++k) {
554 qij += X_(k,i) * X_(k,j);
555 }
556 err = (i==j ? std::abs(qij-one) : std::abs(qij));
557 maxerr = (err > maxerr ? err : maxerr);
558 if (verbosity > 1) {
559 outStream << std::setw(12) << std::right << qij;
560 }
561 }
562 if (verbosity > 1) {
563 outStream << std::endl;
564 }
565 if (maxerr > tol) {
566 break;
567 }
568 }
569 if (verbosity > 0) {
570 outStream << std::endl << " TEST ORTHOGONALIZATION: Max Error = "
571 << std::setw(12) << std::right << maxerr
572 << std::endl;
573 outStream.flags(oflags);
574 }
575 return (maxerr < tol ? true : false);
576 }
577
578}; // class Sketch
579
580} // namespace ROL
581
582#endif
Vector< Real > V
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Contains definitions of custom data types in ROL.
Provides an interface for randomized sketching.
int computeC(void)
LA::Matrix< Real > C_
Ptr< std::mt19937_64 > gen_
std::vector< Ptr< Vector< Real > > > Upsilon_
void mgs2(std::vector< Ptr< Vector< Real > > > &Y) const
LA::Matrix< Real > X_
LA::Matrix< Real > Psi_
virtual ~Sketch(void)
int reconstruct(Vector< Real > &a, const int col)
bool testP(std::ostream &outStream=std::cout, const int verbosity=0)
LA::Matrix< Real > Omega_
int advance(const Real nu, Vector< Real > &h, const int col, const Real eta=1.0)
const Real orthTol_
LAPACK< int, Real > lapack_
Ptr< std::ostream > out_
void update(void)
std::vector< Ptr< Vector< Real > > > Y_
Ptr< std::normal_distribution< Real > > dist_
LA::Matrix< Real > Z_
int computeP(void)
void setStream(Ptr< std::ostream > &out)
int LSsolver(LA::Matrix< Real > &A, LA::Matrix< Real > &B, const bool trans=false) const
Ptr< Elementwise::NormalRandom< Real > > nrand_
const bool truncate_
std::vector< Ptr< Vector< Real > > > Phi_
bool testQ(std::ostream &outStream=std::cout, const int verbosity=0)
int lowRankApprox(LA::Matrix< Real > &A, const int r) const
void setRank(const int rank)
bool test(const int rank, std::ostream &outStream=std::cout, const int verbosity=0)
const int orthIt_
int computeQ(void)
void reset(void)
Sketch(const Vector< Real > &x, const int ncol, const int rank, const Real orthTol=1e-8, const int orthIt=2, const bool truncate=false, const unsigned dom_seed=0, const unsigned rng_seed=0)
Defines the linear algebra or vector space interface.
virtual Real norm() const =0
Returns where .
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
virtual Real dot(const Vector &x) const =0
Compute where .