1 #include <iostream>
2 ///
3 /// With header __flens.cxx__ all of FLENS gets included.
4 ///
5 /// :links: __flens.cxx__ -> file:flens/flens.cxx
6 #include <flens/flens.cxx>
7
8 using namespace std;
9 using namespace flens;
10
11 typedef double T;
12
13 int
14 main()
15 {
16 ///
17 /// Define some convenient typedefs for the matrix/vector types
18 ///
19 typedef GeMatrix<FullStorage<T> > GeMatrix;
20 typedef DenseVector<Array<T> > DenseVector;
21
22 ///
23 /// Define an underscore operator for convenient matrix slicing
24 ///
25 typedef GeMatrix::IndexType IndexType;
26 const Underscore<IndexType> _;
27
28 const IndexType m = 4,
29 n = 5;
30
31 GeMatrix Ab(m, n);
32
33 Ab = 2, 3, -1, 0, 20,
34 -6, -5, 0, 2, -33,
35 2, -5, 6, -6, -43,
36 4, 6, 2, -3, 49;
37
38 cout << "Ab = " << Ab << endl;
39
40 ///
41 /// `tau` will contain the scalar factors of the elementary reflectors
42 /// (see __dgeqrf__ for details). Vector `work` is used as workspace
43 /// and, if empty, gets resized by __lapack::qrf__ to optimal size.
44 ///
45 DenseVector tau(min(m,n));
46 DenseVector work;
47
48 ///
49 /// Compute the $QR$ factorization of $A$ with __lapack::qrf__ the
50 /// FLENS port of LAPACK's __dgeqrf__.
51 ///
52 lapack::qrf(Ab, tau, work);
53 cout << "Ab = " << Ab << endl;
54 cout << "tau = " << tau << endl;
55
56 ///
57 /// Solve the system of linear equation $Ax =B$ using the triangular
58 /// solver __blas::sm__ which is the FLENS implementation of the BLAS
59 /// routine __trsm__. Note that `A.upper()` creates a triangular view.
60 ///
61 const auto A = Ab(_,_(1,m));
62 auto B = Ab(_,_(m+1,n));
63
64 blas::sm(Left, NoTrans, 1, A.upper(), B);
65
66 cout << "X = " << B << endl;
67 }
68
69 ///
70 /// :links: __lapack::qrf__ -> file:flens/lapack/qr/qrf.h
71 /// __dgeqrf__ -> file:flens/lapack/interface/ref_lapack/dgeqrf.f
72 /// __blas::sm__ -> file:flens/blas/level3/sm.h
73 /// __trsm__ -> file:flens/lapack/interface/ref_blas/dtrsm.f
74 ///
2 ///
3 /// With header __flens.cxx__ all of FLENS gets included.
4 ///
5 /// :links: __flens.cxx__ -> file:flens/flens.cxx
6 #include <flens/flens.cxx>
7
8 using namespace std;
9 using namespace flens;
10
11 typedef double T;
12
13 int
14 main()
15 {
16 ///
17 /// Define some convenient typedefs for the matrix/vector types
18 ///
19 typedef GeMatrix<FullStorage<T> > GeMatrix;
20 typedef DenseVector<Array<T> > DenseVector;
21
22 ///
23 /// Define an underscore operator for convenient matrix slicing
24 ///
25 typedef GeMatrix::IndexType IndexType;
26 const Underscore<IndexType> _;
27
28 const IndexType m = 4,
29 n = 5;
30
31 GeMatrix Ab(m, n);
32
33 Ab = 2, 3, -1, 0, 20,
34 -6, -5, 0, 2, -33,
35 2, -5, 6, -6, -43,
36 4, 6, 2, -3, 49;
37
38 cout << "Ab = " << Ab << endl;
39
40 ///
41 /// `tau` will contain the scalar factors of the elementary reflectors
42 /// (see __dgeqrf__ for details). Vector `work` is used as workspace
43 /// and, if empty, gets resized by __lapack::qrf__ to optimal size.
44 ///
45 DenseVector tau(min(m,n));
46 DenseVector work;
47
48 ///
49 /// Compute the $QR$ factorization of $A$ with __lapack::qrf__ the
50 /// FLENS port of LAPACK's __dgeqrf__.
51 ///
52 lapack::qrf(Ab, tau, work);
53 cout << "Ab = " << Ab << endl;
54 cout << "tau = " << tau << endl;
55
56 ///
57 /// Solve the system of linear equation $Ax =B$ using the triangular
58 /// solver __blas::sm__ which is the FLENS implementation of the BLAS
59 /// routine __trsm__. Note that `A.upper()` creates a triangular view.
60 ///
61 const auto A = Ab(_,_(1,m));
62 auto B = Ab(_,_(m+1,n));
63
64 blas::sm(Left, NoTrans, 1, A.upper(), B);
65
66 cout << "X = " << B << endl;
67 }
68
69 ///
70 /// :links: __lapack::qrf__ -> file:flens/lapack/qr/qrf.h
71 /// __dgeqrf__ -> file:flens/lapack/interface/ref_lapack/dgeqrf.f
72 /// __blas::sm__ -> file:flens/blas/level3/sm.h
73 /// __trsm__ -> file:flens/lapack/interface/ref_blas/dtrsm.f
74 ///