Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
numerics/example/DenseMatrix/cxx_main.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Teuchos: Common Tools Package
6// Copyright (2004) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41// @HEADER
42*/
43
47#include "Teuchos_RCP.hpp"
48#include "Teuchos_Version.hpp"
49
50int main(int argc, char* argv[])
51{
52 std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
53
54 // Creating a double-precision matrix can be done in several ways:
55 // Create an empty matrix with no dimension
57 // Create an empty 3x4 matrix
59 // Basic copy of My_Matrix
60 Teuchos::SerialDenseMatrix<int,double> My_Copy1( My_Matrix ),
61 // (Deep) Copy of principle 3x3 submatrix of My_Matrix
62 My_Copy2( Teuchos::Copy, My_Matrix, 3, 3 ),
63 // (Shallow) Copy of 2x3 submatrix of My_Matrix
64 My_Copy3( Teuchos::View, My_Matrix, 2, 3, 1, 1 );
65 // Create a double-precision vector:
67
68 // The matrix dimensions and strided storage information can be obtained:
69 int rows, cols, stride;
70 rows = My_Copy3.numRows(); // number of rows
71 cols = My_Copy3.numCols(); // number of columns
72 stride = My_Copy3.stride(); // storage stride
75 TEUCHOS_ASSERT_EQUALITY(stride, 3);
76
77 // Matrices can change dimension:
78 Empty_Matrix.shape( 3, 3 ); // size non-dimensional matrices
79 My_Matrix.reshape( 3, 3 ); // resize matrices and save values
80
81 // Filling matrices with numbers can be done in several ways:
82 My_Matrix.random(); // random numbers
83 My_Copy1.putScalar( 1.0 ); // every entry is 1.0
84 My_Copy2(1,1) = 10.0; // individual element access
85 Empty_Matrix = My_Matrix; // copy My_Matrix to Empty_Matrix
86 x = 1.0; // every entry of vector is 1.0
87 y = 1.0;
88
89 // Basic matrix arithmetic can be performed:
90 double d;
92 // Matrix multiplication ( My_Prod = 1.0*My_Matrix*My_Copy^T )
94 1.0, My_Matrix, My_Copy3, 0.0 );
95 My_Copy2 += My_Matrix; // Matrix addition
96 My_Copy2.scale( 0.5 ); // Matrix scaling
97 d = x.dot( y ); // Vector dot product
98 (void)d; // Not used!
99
100 // The pointer to the array of matrix values can be obtained:
101 double *My_Array=0, *My_Column=0;
102 My_Array = My_Matrix.values(); // pointer to matrix values
103 My_Column = My_Matrix[2]; // pointer to third column values
104 (void)My_Array; // Not used!
105 (void)My_Column; // Not used!
106
107 // The norm of a matrix can be computed:
108 double norm_one, norm_inf, norm_fro;
109 norm_one = My_Matrix.normOne(); // one norm
110 norm_inf = My_Matrix.normInf(); // infinity norm
111 norm_fro = My_Matrix.normFrobenius(); // frobenius norm
112 (void)norm_one; // Not used!
113 (void)norm_inf; // Not used!
114 (void)norm_fro; // Not used!
115
116 // Matrices can be compared:
117 // Check if the matrices are equal in dimension and values
118 if (Empty_Matrix == My_Matrix) {
119 std::cout<< "The matrices are the same!" <<std::endl;
120 }
121 // Check if the matrices are different in dimension or values
122 if (My_Copy2 != My_Matrix) {
123 std::cout<< "The matrices are different!" <<std::endl;
124 }
125
126 // A matrix can be factored and solved using Teuchos::SerialDenseSolver.
129 X.putScalar(1.0);
130 B.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, My_Matrix, X, 0.0 );
131 X.putScalar(0.0); // Make sure the computed answer is correct.
132
133 int info = 0;
134 My_Solver.setMatrix( Teuchos::rcp( &My_Matrix, false ) );
135 My_Solver.setVectors( Teuchos::rcp( &X, false ), Teuchos::rcp( &B, false ) );
136 info = My_Solver.factor();
137 if (info != 0)
138 std::cout << "Teuchos::SerialDenseSolver::factor() returned : " << info << std::endl;
139 info = My_Solver.solve();
140 if (info != 0)
141 std::cout << "Teuchos::SerialDenseSolver::solve() returned : " << info << std::endl;
142
143 // A matrix can be sent to the output stream:
144 std::cout<< std::endl << printMat(My_Matrix) << std::endl;
145 std::cout<< printMat(X) << std::endl;
146
147 return 0;
148}
Reference-counted pointer class and non-member templated function implementations.
Templated serial dense matrix class.
Templated class for solving dense linear problems.
Templated serial dense vector class.
This class creates and provides basic support for dense rectangular matrix of templated type.
ScalarTraits< ScalarType >::magnitudeType normOne() const
Returns the 1-norm of the matrix.
int reshape(OrdinalType numRows, OrdinalType numCols)
Reshaping method for changing the size of a SerialDenseMatrix, keeping the entries.
int scale(const ScalarType alpha)
Scale this matrix by alpha; *this = alpha**this.
OrdinalType stride() const
Returns the stride between the columns of this matrix in memory.
int random()
Set all values in the matrix to be random numbers.
ScalarTraits< ScalarType >::magnitudeType normInf() const
Returns the Infinity-norm of the matrix.
OrdinalType numRows() const
Returns the row dimension of this matrix.
ScalarType * values() const
Data array access method.
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Multiply A * B and add them to this; this = beta * this + alpha*A*B.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
OrdinalType numCols() const
Returns the column dimension of this matrix.
int shape(OrdinalType numRows, OrdinalType numCols)
Shape method for changing the size of a SerialDenseMatrix, initializing entries to zero.
A class for solving dense linear problems.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
int factor()
Computes the in-place LU factorization of the matrix using the LAPACK routine _GETRF.
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
This class creates and provides basic support for dense vectors of templated type as a specialization...
int main()
Definition: evilMain.cpp:75
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
This macro is checks that to numbers are equal and if not then throws an exception with a good error ...
Definition: PackageB.cpp:3
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
std::string Teuchos_Version()