本文整理汇总了C++中teuchos::LAPACK::TRTRS方法的典型用法代码示例。如果您正苦于以下问题：C++ LAPACK::TRTRS方法的具体用法？C++ LAPACK::TRTRS怎么用？C++ LAPACK::TRTRS使用的例子？那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类teuchos::LAPACK的用法示例。
在下文中一共展示了LAPACK::TRTRS方法的3个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: gmres
//.........这里部分代码省略.........
//set vi to be ith col of V
Teuchos::SerialDenseMatrix<int, double> vi(Teuchos::Copy, V, n,1,0,i);
//Calculate inner product
ip.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, vi, w, 0.0);
h(i,k-1)= ip(0,0);
//scale vi by h(i,k-1)
vi.scale(ip(0,0));
w-=vi;
}
h(k,k-1)=w.normFrobenius();
w.scale(1.0/w.normFrobenius());
//add column vk+1=w to V
for (int i=0; i<n; i++){
V(i,k)=w(i,0);
}
//Solve upper hessenberg least squares problem via Givens rotations
//Compute previous Givens rotations
for (int i=0; i<k-1; i++){
// double hi=h(i,k-1);
// double hi1=h(i+1,k-1);
// h(i,k-1)=c(i,0)*h(i,k-1)+s(i,0)*h(i+1,k-1);
// h(i+1,k-1)=-1*s(i,0)*h(i,k-1)+c(i,0)*h(i+1,k-1);
// h(i,k-1)=c(i,0)*hi+s(i,0)*hi1;
// h(i+1,k-1)=-1*s(i,0)*hi+c(i,0)*hi1;
double q=c(i,0)*h(i,k-1)+s(i,0)*h(i+1,k-1);
h(i+1,k-1)=-1*s(i,0)*h(i,k-1)+c(i,0)*h(i+1,k-1);
h(i,k-1)=q;
}
//Compute next Givens rotations
c.reshape(k,1);
s.reshape(k,1);
bb.reshape(k+1,1);
double l = sqrt(h(k-1,k-1)*h(k-1,k-1)+h(k,k-1)*h(k,k-1));
c(k-1,0)=h(k-1,k-1)/l;
s(k-1,0)=h(k,k-1)/l;
std::cout << "c " << c(k-1,0)<<std::endl;
std::cout << "s " << s(k-1,0)<<std::endl;
// Givens rotation on h and bb
// h(k-1,k-1)=l;
// h(k,k-1)=0;
double hk=h(k,k-1);
double hk1=h(k-1,k-1);
h(k-1,k-1)=c(k-1,0)*hk1+s(k-1,0)*hk;
h(k,k-1)=-1*s(k-1,0)*hk1+c(k-1,0)*hk;
std::cout << "l = " << l <<std::endl;
std::cout << "h(k-1,k-1) = should be l " << h(k-1,k-1) <<std::endl;
std::cout << "h(k,k-1) = should be 0 " << h(k,k-1) <<std::endl;
bb(k,0)=-1*s(k-1,0)*bb(k-1,0);
bb(k-1,0)=c(k-1,0)*bb(k-1,0);
//Determine residual
resid =fabs(bb(k,0));
std::cout << "resid = " << resid <<std::endl;
k++;
}
//Extract upper triangular square matrix
bb.reshape(h.numRows()-1 ,1);
//Solve linear system
int info;
std::cout << "bb pre solve = " << bb << std::endl;
std::cout << "h= " << h << std::endl;
Teuchos::LAPACK<int, double> lapack;
lapack.TRTRS('U', 'N', 'N', h.numRows()-1, 1, h.values(), h.stride(), bb.values(), bb.stride(),&info);
V.reshape(n,k-1);
std::cout << "V= " << V << std::endl;
std::cout << "y= " << bb << std::endl;
X.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, V, bb, 1.0);
std::cout << "X= " << X << std::endl;
//Check V is orthogoanl
// Teuchos::SerialDenseMatrix<int, double> vtv(V);
// vtv.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, V, V, 0.0);
// std::cout << "Vtv" << vtv << std::endl;
return 0;
}示例2: precond
//.........这里部分代码省略.........
V.reshape(n,k+1);
//set vk to be kth col of V
Teuchos::SerialDenseMatrix<int, double> vk(Teuchos::Copy, V, n,1,0,k-1);
//Preconditioning step: solve Mz=vk
Teuchos::SerialDenseMatrix<int, double> z(vk);
if (PrecNum == 1) {
Stokhos::DiagPreconditioner precond(M);
precond.ApplyInverse(vk,z,prec_iter);
}
else if (PrecNum == 2) {
Stokhos::JacobiPreconditioner precond(M);
precond.ApplyInverse(vk,z,2);
}
else if (PrecNum == 3) {
Stokhos::GSPreconditioner precond(M,1);
precond.ApplyInverse(vk,z,1);
}
else if (PrecNum == 4) {
Stokhos::SchurPreconditioner precond(M, order, dim, diag);
precond.ApplyInverse(vk,z,prec_iter);
}
w.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1, A, z, 0.0);
Teuchos::SerialDenseMatrix<int, double> vi(n,1);
Teuchos::SerialDenseMatrix<int, double> ip(1,1);
for (int i=0; i<k; i++) {
//set vi to be ith col of V
Teuchos::SerialDenseMatrix<int, double> vi(Teuchos::Copy, V, n,1,0,i);
//Calculate inner product
ip.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, vi, w, 0.0);
h(i,k-1)= ip(0,0);
//scale vi by h(i,k-1)
vi.scale(ip(0,0));
w-=vi;
}
h(k,k-1)=w.normFrobenius();
w.scale(1.0/h(k,k-1));
//add column vk+1=w to V
for (int i=0; i<n; i++) {
V(i,k)=w(i,0);
}
//Solve upper hessenberg least squares problem via Givens rotations
//Compute previous Givens rotations
for (int i=0; i<k-1; i++) {
double q=c(i,0)*h(i,k-1)+s(i,0)*h(i+1,k-1);
h(i+1,k-1)=-1*s(i,0)*h(i,k-1)+c(i,0)*h(i+1,k-1);
h(i,k-1)=q;
}
//Compute next Givens rotations
c.reshape(k,1);
s.reshape(k,1);
bb.reshape(k+1,1);
double l = sqrt(h(k-1,k-1)*h(k-1,k-1)+h(k,k-1)*h(k,k-1));
c(k-1,0)=h(k-1,k-1)/l;
s(k-1,0)=h(k,k-1)/l;
// Givens rotation on h and bb
h(k-1,k-1)=l;
h(k,k-1)=0;
bb(k,0)=-s(k-1,0)*bb(k-1,0);
bb(k-1,0)=c(k-1,0)*bb(k-1,0);
//Determine residual
resid = fabs(bb(k,0));
k++;
}
//Extract upper triangular square matrix
bb.reshape(h.numRows()-1 ,1);
//Solve linear system
int info;
Teuchos::LAPACK<int, double> lapack;
lapack.TRTRS('U', 'N', 'N', h.numRows()-1, 1, h.values(), h.stride(), bb.values(), bb.stride(),&info);
Teuchos::SerialDenseMatrix<int, double> ans(X);
V.reshape(n,k-1);
ans.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, V, bb, 0.0);
if (PrecNum == 1) {
Stokhos::DiagPreconditioner precond(M);
precond.ApplyInverse(ans,ans,prec_iter);
}
else if (PrecNum == 2) {
Stokhos::JacobiPreconditioner precond(M);
precond.ApplyInverse(ans,ans,2);
}
else if (PrecNum == 3) {
Stokhos::GSPreconditioner precond(M,1);
precond.ApplyInverse(ans,ans,1);
}
else if (PrecNum == 4) {
Stokhos::SchurPreconditioner precond(M, order, dim, diag);
precond.ApplyInverse(ans,ans,prec_iter);
}
X+=ans;
std::cout << "iteration count= " << k-1 << std::endl;
return 0;
}示例3: pregmres
//.........这里部分代码省略.........
r0.scale(1/resid);
Teuchos::SerialDenseMatrix<int, double> h(1,1);
//Matrix of orthog basis vectors V
Teuchos::SerialDenseMatrix<int, double> V(n,1);
//Set v=r0/norm(r0) to be 1st col of V
for (int i=0; i<n; i++){
V(i,0)=r0(i,0);
}
//right hand side
Teuchos::SerialDenseMatrix<int, double> bb(1,1);
bb(0,0)=resid;
Teuchos::SerialDenseMatrix<int, double> w(n,1);
Teuchos::SerialDenseMatrix<int, double> c;
Teuchos::SerialDenseMatrix<int, double> s;
while (resid > tolerance && k < max_iter){
std::cout << "k = " << k << std::endl;
h.reshape(k+1,k);
//Arnoldi iteration(Gram-Schmidt )
V.reshape(n,k+1);
//set vk to be kth col of V
Teuchos::SerialDenseMatrix<int, double> vk(Teuchos::Copy, V, n,1,0,k-1);
//Preconditioning step w=AMj(-1)vj
w.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1/D(k-1,0), A, vk, 0.0);
Teuchos::SerialDenseMatrix<int, double> vi(n,1);
Teuchos::SerialDenseMatrix<int, double> ip(1,1);
for (int i=0; i<k; i++){
//set vi to be ith col of V
Teuchos::SerialDenseMatrix<int, double> vi(Teuchos::Copy, V, n,1,0,i);
//Calculate inner product
ip.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, vi, w, 0.0);
h(i,k-1)= ip(0,0);
//scale vi by h(i,k-1)
vi.scale(ip(0,0));
w-=vi;
}
h(k,k-1)=w.normFrobenius();
w.scale(1.0/w.normFrobenius());
//add column vk+1=w to V
for (int i=0; i<n; i++){
V(i,k)=w(i,0);
}
//Solve upper hessenberg least squares problem via Givens rotations
//Compute previous Givens rotations
for (int i=0; i<k-1; i++){
h(i,k-1)=c(i,0)*h(i,k-1)+s(i,0)*h(i+1,k-1);
h(i+1,k-1)=-s(i,0)*h(i,k-1)+c(i,0)*h(i+1,k-1);
}
//Compute next Givens rotations
c.reshape(k,1);
s.reshape(k,1);
bb.reshape(k+1,1);
double l = sqrt(h(k-1,k-1)*h(k-1,k-1)+h(k,k-1)*h(k,k-1));
c(k-1,0)=h(k-1,k-1)/l;
s(k-1,0)=h(k,k-1)/l;
std::cout <<" h(k,k-1)= " << h(k,k-1) << std::endl;
// Givens rotation on h and bb
h(k-1,k-1)=l;
h(k,k-1)=0;
bb(k-1,0)=c(k-1,0)*bb(k-1,0);
bb(k,0)=-s(k-1,0)*bb(k-1,0);
//Determine residual
resid = fabs(bb(k,0));
std::cout << "resid = " << resid << std::endl;
k=k+1;
}
//Extract upper triangular square matrix
bb.reshape(h.numRows()-1 ,1);
//Solve linear system
int info;
Teuchos::LAPACK<int, double> lapack;
lapack.TRTRS('U', 'N', 'N', h.numRows()-1, 1, h.values(), h.stride(), bb.values(), bb.stride(),&info);
//Found y=Mx
for (int i=0; i<k-1; i++){
bb(i,0)=bb(i,0)/D(i,0);
}
V.reshape(n,k-1);
Teuchos::SerialDenseMatrix<int, double> ans(X);
ans.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, V, bb, 1.0);
std::cout << "ans= " << ans << std::endl;
std::cout << "h= " << h << std::endl;
return 0;
}