本文整理汇总了C++中UpperTriangularMatrix类的典型用法代码示例。如果您正苦于以下问题:C++ UpperTriangularMatrix类的具体用法?C++ UpperTriangularMatrix怎么用?C++ UpperTriangularMatrix使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了UpperTriangularMatrix类的14个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: QRZ
void QRZ(Matrix& X, UpperTriangularMatrix& U)
{
REPORT
Tracer et("QRZ(1)");
int n = X.Nrows(); int s = X.Ncols(); U.ReSize(s); U = 0.0;
Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
int j, k; int J = s; int i = s;
while (i--)
{
Real* xj0 = xi0; Real* xi = xi0; k = n;
if (k) for (;;)
{
u = u0; Real Xi = *xi; Real* xj = xj0;
j = J; while(j--) *u++ += Xi * *xj++;
if (!(--k)) break;
xi += s; xj0 += s;
}
Real sum = sqrt(*u0); *u0 = sum; u = u0+1;
if (sum == 0.0) Throw(SingularException(U));
int J1 = J-1; j = J1; while(j--) *u++ /= sum;
xj0 = xi0; xi = xi0++; k = n;
if (k) for (;;)
{
u = u0+1; Real Xi = *xi; Real* xj = xj0;
Xi /= sum; *xj++ = Xi;
j = J1; while(j--) *xj++ -= *u++ * Xi;
if (!(--k)) break;
xi += s; xj0 += s;
}
u0 += J--;
}
}
示例2: MakeCovariance
void NonLinearLeastSquares::MakeCovariance()
{
if (Covariance.Nrows()==0)
{
UpperTriangularMatrix UI = U.i();
Covariance << UI * UI.t() * errorvar;
SE << Covariance; // get diagonals
for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
}
}
示例3: M
inline
Matrix<C>::Matrix(const UpperTriangularMatrix<C>& M)
: entries_(M.row_dimension()*M.column_dimension()),
rowdim_(M.row_dimension()), coldim_(M.column_dimension())
{
for (size_type i(0); i < rowdim_; i++)
for (size_type j(i); j < coldim_; j++)
{
this->operator () (i, j) = M(i, j);
}
}
示例4: updateQRZ
void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U)
{
REPORT
Tracer et("updateQRZ(3)");
int s = X.Ncols();
if (s != U.Ncols())
Throw(ProgramException("Incompatible dimensions",X,U));
if (s == 0) return;
Real* xi0 = X.data(); Real* u = U.data();
for (int i=1; i<=s; ++i)
{
Real r = *u; Real sum = 0.0;
{
Real* xi=xi0; int k=i; int l=s;
while(k--) { sum += square(*xi); xi+= --l;}
}
sum = sqrt(sum + square(r));
if (sum == 0.0) { REPORT X.column(i) = 0.0; *u = 0.0; }
else
{
Real frs = fabs(r) + sum;
Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
if (r <= 0) { REPORT *u = sum; alpha = -alpha; }
else { REPORT *u = -sum; }
{
Real* xj0=xi0; int k=i; int l=s;
while(k--) { *xj0 *= alpha; --l; xj0 += l;}
}
Real* xj0=xi0; Real* uj=u;
for (int j=i+1; j<=s; ++j)
{
Real sum = 0.0; ++xj0; ++uj;
Real* xi=xi0; Real* xj=xj0; int k=i; int l=s;
while(k--) { sum += *xi * *xj; --l; xi += l; xj += l; }
sum += a0 * *uj;
xi=xi0; xj=xj0; k=i; l=s;
while(k--) { *xj -= sum * *xi; --l; xi += l; xj += l; }
*uj -= sum * a0;
}
}
++xi0; u += s-i+1;
}
}
示例5: downdate_Cholesky
// produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x)
{
int nRC = chol.Nrows();
// solve R^T a = x
LowerTriangularMatrix L = chol.t();
ColumnVector a(nRC); a = 0.0;
int i, j;
for (i = 1; i <= nRC; ++i)
{
// accumulate subtr sum
Real subtrsum = 0.0;
for(int k = 1; k < i; ++k) subtrsum += a(k) * L(i,k);
a(i) = (x(i) - subtrsum) / L(i,i);
}
// test that l2 norm of a is < 1
Real squareNormA = a.SumSquare();
if (squareNormA >= 1.0)
Throw(ProgramException("downdate_Cholesky() fails", chol));
Real alpha = sqrt(1.0 - squareNormA);
// compute and apply Givens rotations to the vector a
ColumnVector cGivens(nRC); cGivens = 0.0;
ColumnVector sGivens(nRC); sGivens = 0.0;
for(i = nRC; i >= 1; i--)
alpha = pythag(alpha, a(i), cGivens(i), sGivens(i));
// apply Givens rotations to the jth column of chol
ColumnVector xtilde(nRC); xtilde = 0.0;
for(j = nRC; j >= 1; j--)
{
// only the first j rotations have an affect on chol,0
for(int k = j; k >= 1; k--)
GivensRotation(cGivens(k), -sGivens(k), chol(k,j), xtilde(j));
}
}
示例6: left_circular_update_Cholesky
// produces the Cholesky decomposition of EAE where A = chol.t() * chol
// and E produces a LEFT circular shift of the rows and columns from
// 1,...,k-1,k,k+1,...l,l+1,...,p to
// 1,...,k-1,k+1,...l,k,l+1,...,p to
void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
{
int nRC = chol.Nrows();
int i, j;
// I. compute shift of column k to the lth position
Matrix cholCopy = chol;
// a. grab column k
ColumnVector columnK = cholCopy.Column(k);
// b. shift columns k+1,...l to the LEFT
for(j = k+1; j <= l; ++j)
cholCopy.Column(j-1) = cholCopy.Column(j);
// c. copy the elements of columnK into the lth column of cholCopy
cholCopy.Column(l) = 0.0;
for(i = 1; i <= k; ++i)
cholCopy(i,l) = columnK(i);
// II. apply and compute Given's rotations
int nGivens = l-k;
ColumnVector cGivens(nGivens); cGivens = 0.0;
ColumnVector sGivens(nGivens); sGivens = 0.0;
for(j = k; j <= nRC; ++j)
{
ColumnVector columnJ = cholCopy.Column(j);
// apply the previous Givens rotations to columnJ
int imax = j - k; if (imax > nGivens) imax = nGivens;
for(int i = 1; i <= imax; ++i)
{
int gIndex = i;
int topRowIndex = k + i - 1;
GivensRotationR(cGivens(gIndex), sGivens(gIndex),
columnJ(topRowIndex), columnJ(topRowIndex+1));
}
// compute a new Given's rotation when j < l
if(j < l)
{
int gIndex = j-k+1;
columnJ(j) = pythag(columnJ(j), columnJ(j+1), cGivens(gIndex),
sGivens(gIndex));
columnJ(j+1) = 0.0;
}
cholCopy.Column(j) = columnJ;
}
chol << cholCopy;
}
示例7: test4
void test4(Real* y, Real* x1, Real* x2, int nobs, int npred)
{
cout << "\n\nTest 4 - QR triangularisation\n";
// QR triangularisation method
// load data - 1s into col 1 of matrix
int npred1 = npred+1;
Matrix X(nobs,npred1); ColumnVector Y(nobs);
X.Column(1) = 1.0; X.Column(2) << x1; X.Column(3) << x2; Y << y;
// do Householder triangularisation
// no need to deal with constant term separately
Matrix X1 = X; // Want copy of matrix
ColumnVector Y1 = Y;
UpperTriangularMatrix U; ColumnVector M;
QRZ(X1, U); QRZ(X1, Y1, M); // Y1 now contains resids
ColumnVector A = U.i() * M;
ColumnVector Fitted = X * A;
Real ResVar = Y1.SumSquare() / (nobs-npred1);
// get variances of estimates
U = U.i(); DiagonalMatrix D; D << U * U.t();
// Get diagonals of Hat matrix
DiagonalMatrix Hat; Hat << X1 * X1.t();
// print out answers
cout << "\nEstimates and their standard errors\n\n";
ColumnVector SE(npred1);
for (int i=1; i<=npred1; i++) SE(i) = sqrt(D(i)*ResVar);
cout << setw(11) << setprecision(5) << (A | SE) << endl;
cout << "\nObservations, fitted value, residual value, hat value\n";
cout << setw(9) << setprecision(3) <<
(X.Columns(2,3) | Y | Fitted | Y1 | Hat.AsColumn());
cout << "\n\n";
}
示例8: right_circular_update_Cholesky
// produces the Cholesky decomposition of EAE where A = chol.t() * chol
// and E produces a RIGHT circular shift of the rows and columns from
// 1,...,k-1,k,k+1,...l,l+1,...,p to
// 1,...,k-1,l,k,k+1,...l-1,l+1,...p
void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
{
int nRC = chol.Nrows();
int i, j;
// I. compute shift of column l to the kth position
Matrix cholCopy = chol;
// a. grab column l
ColumnVector columnL = cholCopy.Column(l);
// b. shift columns k,...l-1 to the RIGHT
for(j = l-1; j >= k; --j)
cholCopy.Column(j+1) = cholCopy.Column(j);
// c. copy the top k-1 elements of columnL into the kth column of cholCopy
cholCopy.Column(k) = 0.0;
for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);
// II. determine the l-k Given's rotations
int nGivens = l-k;
ColumnVector cGivens(nGivens); cGivens = 0.0;
ColumnVector sGivens(nGivens); sGivens = 0.0;
for(i = l; i > k; i--)
{
int givensIndex = l-i+1;
columnL(i-1) = pythag(columnL(i-1), columnL(i),
cGivens(givensIndex), sGivens(givensIndex));
columnL(i) = 0.0;
}
// the kth entry of columnL is the new diagonal element in column k of cholCopy
cholCopy(k,k) = columnL(k);
// III. apply these Given's rotations to subsequent columns
// for columns k+1,...,l-1 we only need to apply the last nGivens-(j-k) rotations
for(j = k+1; j <= nRC; ++j)
{
ColumnVector columnJ = cholCopy.Column(j);
int imin = nGivens - (j-k) + 1; if (imin < 1) imin = 1;
for(int gIndex = imin; gIndex <= nGivens; ++gIndex)
{
// apply gIndex Given's rotation
int topRowIndex = k + nGivens - gIndex;
GivensRotationR(cGivens(gIndex), sGivens(gIndex),
columnJ(topRowIndex), columnJ(topRowIndex+1));
}
cholCopy.Column(j) = columnJ;
}
chol << cholCopy;
}
示例9: Print
void Print(const UpperTriangularMatrix& X)
{
++PCN;
cout << "\nMatrix type: " << X.Type().Value() << " (";
cout << X.Nrows() << ", ";
cout << X.Ncols() << ")\n\n";
if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
int nr=X.Nrows(); int nc=X.Ncols();
for (int i=1; i<=nr; i++)
{
int j;
for (j=1; j<i; j++) cout << "\t";
for (j=i; j<=nc; j++) cout << X(i,j) << "\t";
cout << "\n";
}
cout << flush; ++PCZ;
}
示例10: update_Cholesky
// produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
void update_Cholesky(UpperTriangularMatrix &chol, RowVector x)
{
int nc = chol.Nrows();
ColumnVector cGivens(nc); cGivens = 0.0;
ColumnVector sGivens(nc); sGivens = 0.0;
for(int j = 1; j <= nc; ++j) // process the jth column of chol
{
// apply the previous Givens rotations k = 1,...,j-1 to column j
for(int k = 1; k < j; ++k)
GivensRotation(cGivens(k), sGivens(k), chol(k,j), x(j));
// determine the jth Given's rotation
pythag(chol(j,j), x(j), cGivens(j), sGivens(j));
// apply the jth Given's rotation
{
Real tmp0 = cGivens(j) * chol(j,j) + sGivens(j) * x(j);
chol(j,j) = tmp0; x(j) = 0.0;
}
}
}
示例11: updateQRZ
void updateQRZ(Matrix& X, UpperTriangularMatrix& U)
{
REPORT
Tracer et("updateQRZ");
int n = X.Nrows(); int s = X.Ncols();
if (s != U.Ncols())
Throw(ProgramException("Incompatible dimensions",X,U));
if (n == 0 || s == 0) return;
Real* xi0 = X.Store(); Real* u0 = U.Store(); Real* u;
RowVector V(s); Real* v0 = V.Store(); Real* v; V = 0.0;
int j, k; int J = s; int i = s;
while (i--)
{
Real* xj0 = xi0; Real* xi = xi0; k = n;
if (k) for (;;)
{
v = v0; Real Xi = *xi; Real* xj = xj0;
j = J; while(j--) *v++ += Xi * *xj++;
if (!(--k)) break;
xi += s; xj0 += s;
}
Real r = *u0;
Real sum = sqrt(*v0 + square(r));
if (sum == 0.0)
{
REPORT
u = u0; v = v0;
j = J; while(j--) { *u++ = 0.0; *v++ = 0.0; }
xj0 = xi0++; k = n;
if (k) for (;;)
{
*xj0 = 0.0;
if (!(--k)) break;
xj0 += s;
}
u0 += J--;
}
else
{
Real frs = fabs(r) + sum;
Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
if (r <= 0) { REPORT alpha = -alpha; *u0 = sum; }
else { REPORT *u0 = -sum; }
j = J - 1; v = v0 + 1; u = u0 + 1;
while (j--)
{ *v = a0 * *u + alpha * *v; *u -= a0 * *v; ++v; ++u; }
xj0 = xi0; xi = xi0++; k = n;
if (k) for (;;)
{
v = v0 + 1; Real Xi = *xi; Real* xj = xj0;
Xi *= alpha; *xj++ = Xi;
j = J - 1; while(j--) *xj++ -= *v++ * Xi;
if (!(--k)) break;
xi += s; xj0 += s;
}
j = J; v = v0;
while (j--) *v++ = 0.0;
u0 += J--;
}
}
}
示例12: trymatc
void trymatc()
{
// cout << "\nTwelfth test of Matrix package\n";
Tracer et("Twelfth test of Matrix package");
Tracer::PrintTrace();
DiagonalMatrix D(15); D=1.5;
Matrix A(15,15);
int i,j;
for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
{ A = A + D; }
ColumnVector B(15);
for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
{
Tracer et1("Stage 1");
ColumnVector B1=B;
B=(A*2.0).i() * B1;
Matrix X = A*B-B1/2.0;
Clean(X, 0.000000001); Print(X);
A.ReSize(3,5);
for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
B = A.AsColumn()+10000;
RowVector R = (A+10000).AsColumn().t();
Print( RowVector(R-B.t()) );
}
{
Tracer et1("Stage 2");
B = A.AsColumn()+10000;
Matrix XR = (A+10000).AsMatrix(15,1).t();
Print( RowVector(XR-B.t()) );
}
{
Tracer et1("Stage 3");
B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
Matrix MR = (A+10000).AsColumn().t();
Print( RowVector(MR-B.t()) );
B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
MR = A.AsColumn().t();
Print( RowVector(MR-B.t()) );
}
{
Tracer et1("Stage 4");
B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
RowVector R = A.AsColumn().t();
Print( RowVector(R-B.t()) );
}
{
Tracer et1("Stage 5");
RowVector R = (A.AsColumn()-5000).t();
B = ((R.t()+10000) - A.AsColumn())-5000;
Print( RowVector(B.t()) );
}
{
Tracer et1("Stage 6");
B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
Print(ColumnVector(B1-B));
}
{
Tracer et1("Stage 7");
Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
Print(B);
}
{
Tracer et1("Stage 8");
A.ReSize(7,7); D.ReSize(7);
for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
for (i=1; i<=7; i++) D(i,i) = i;
UpperTriangularMatrix U; U << A;
Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
A.Inject(D); Print(Matrix(X-A));
X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
Print(Matrix(X-A));
}
{
Tracer et1("Stage 9");
A.ReSize(7,5);
for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
Matrix X = A; // X.Release();
Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
}
{
Tracer et1("Stage 10");
// some tests on submatrices
UpperTriangularMatrix U(20);
for (i=1; i<=20; i++) for (j=i; j<=20; j++) U(i,j)=100 * i + j;
UpperTriangularMatrix V = U.SymSubMatrix(1,5);
UpperTriangularMatrix U1 = U;
//.........这里部分代码省略.........
示例13: GEODE_ASSERT
template<class T,int d> void StrainMeasure<T,d>::initialize_rest_state_to_equilateral(const T side_length) {
GEODE_ASSERT(Dm_inverse.size()==elements.size());
UpperTriangularMatrix<T,d> Dm = side_length*equilateral_Dm(Vector<T,d>());
Dm_inverse.fill(Dm.inverse());
}
示例14: trymatc
void trymatc()
{
// cout << "\nTwelfth test of Matrix package\n";
Tracer et("Twelfth test of Matrix package");
Tracer::PrintTrace();
DiagonalMatrix D(15); D=1.5;
Matrix A(15,15);
int i,j;
for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
{ A = A + D; }
ColumnVector B(15);
for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
{
Tracer et1("Stage 1");
ColumnVector B1=B;
B=(A*2.0).i() * B1;
Matrix X = A*B-B1/2.0;
Clean(X, 0.000000001); Print(X);
A.ReSize(3,5);
for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
B = A.AsColumn()+10000;
RowVector R = (A+10000).AsColumn().t();
Print( RowVector(R-B.t()) );
}
{
Tracer et1("Stage 2");
B = A.AsColumn()+10000;
Matrix XR = (A+10000).AsMatrix(15,1).t();
Print( RowVector(XR-B.t()) );
}
{
Tracer et1("Stage 3");
B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
Matrix MR = (A+10000).AsColumn().t();
Print( RowVector(MR-B.t()) );
B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
MR = A.AsColumn().t();
Print( RowVector(MR-B.t()) );
}
{
Tracer et1("Stage 4");
B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
RowVector R = A.AsColumn().t();
Print( RowVector(R-B.t()) );
}
{
Tracer et1("Stage 5");
RowVector R = (A.AsColumn()-5000).t();
B = ((R.t()+10000) - A.AsColumn())-5000;
Print( RowVector(B.t()) );
}
{
Tracer et1("Stage 6");
B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
Print(ColumnVector(B1-B));
}
{
Tracer et1("Stage 7");
Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
Print(B);
}
{
Tracer et1("Stage 8");
A.ReSize(7,7); D.ReSize(7);
for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
for (i=1; i<=7; i++) D(i,i) = i;
UpperTriangularMatrix U; U << A;
Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
A.Inject(D); Print(Matrix(X-A));
X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
Print(Matrix(X-A));
}
{
Tracer et1("Stage 9");
A.ReSize(7,5);
for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
Matrix X = A;
Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
}
{
Tracer et1("Stage 10");
// some tests on submatrices
UpperTriangularMatrix U(20);
for (i=1; i<=20; i++) for (j=i; j<=20; j++) U(i,j)=100 * i + j;
UpperTriangularMatrix V = U.SymSubMatrix(1,5);
UpperTriangularMatrix U1 = U;
//.........这里部分代码省略.........