1 #include <iostream>
2
3 ///
4 /// Include the qd-headers
5 ///
6 #include <qd/qd_real.h>
7 #include <qd/fpu.h>
8
9 #include <flens/flens.cxx>
10
11 using namespace std;
12 using namespace flens;
13
14 ///
15 /// For using quad-double precision in this example change the typedef to
16 ///
17 typedef qd_real T;
18
19 ///
20 /// For understanding the next code lines we first take a look at the
21 /// __QD Library__ documentation:
22 ///
23 /// *--[BOX]------------------------------------------------------------*
24 /// | |
25 /// | The algorithms in the QD library assume IEEE double precision |
26 /// | floating point arithmetic. Since Intel x86 processors have |
27 /// | extended (80-bit) floating point registers, the round-to-double |
28 /// | flag must be enabled in the control word of the FPU for this |
29 /// | library to function properly under x86 processors. The function |
30 /// | `fpu_fix_start` turns on the round-to-double bit in the FPU |
31 /// | control word, while `fpu_fix_end` will restore the original |
32 /// | state. |
33 /// | |
34 /// *-------------------------------------------------------------------*
35 ///
36 /// So the first thing we do in main is turning on the correct rounding ...
37 ///
38 int
39 main()
40 {
41 unsigned int old_cw;
42 fpu_fix_start(&old_cw);
43
44 typedef GeMatrix<FullStorage<T> > Matrix;
45 typedef DenseVector<Array<T> > Vector;
46 typedef Matrix::IndexType IndexType;
47 typedef DenseVector<Array<IndexType> > IndexVector;
48
49 const Underscore<IndexType> _;
50
51 const IndexType m = 4,
52 n = 5;
53
54 Matrix Ab(m, n);
55 IndexVector piv(m);
56
57 Ab = 2, 3, -1, 0, 20,
58 -6, -5, 0, 2, -33,
59 2, -5, 6, -6, -43,
60 4, 6, 2, -3, 49;
61
62 cout << "Ab = " << Ab << endl;
63
64 lapack::trf(Ab, piv);
65
66 const auto A = Ab(_,_(1,m));
67 auto B = Ab(_,_(m+1,n));
68
69 blas::sm(Left, NoTrans, T(1), A.upper(), B);
70
71 cout << "X = " << B << endl;
72
73 ///
74 /// ... and at the end restore FPU rounding behavior as mentioned above.
75 ///
76 fpu_fix_end(&old_cw);
77 }
78
79 ///
80 /// :links: __QD Library__ -> http://crd-legacy.lbl.gov/~dhbailey/mpdist/
81 ///
2
3 ///
4 /// Include the qd-headers
5 ///
6 #include <qd/qd_real.h>
7 #include <qd/fpu.h>
8
9 #include <flens/flens.cxx>
10
11 using namespace std;
12 using namespace flens;
13
14 ///
15 /// For using quad-double precision in this example change the typedef to
16 ///
17 typedef qd_real T;
18
19 ///
20 /// For understanding the next code lines we first take a look at the
21 /// __QD Library__ documentation:
22 ///
23 /// *--[BOX]------------------------------------------------------------*
24 /// | |
25 /// | The algorithms in the QD library assume IEEE double precision |
26 /// | floating point arithmetic. Since Intel x86 processors have |
27 /// | extended (80-bit) floating point registers, the round-to-double |
28 /// | flag must be enabled in the control word of the FPU for this |
29 /// | library to function properly under x86 processors. The function |
30 /// | `fpu_fix_start` turns on the round-to-double bit in the FPU |
31 /// | control word, while `fpu_fix_end` will restore the original |
32 /// | state. |
33 /// | |
34 /// *-------------------------------------------------------------------*
35 ///
36 /// So the first thing we do in main is turning on the correct rounding ...
37 ///
38 int
39 main()
40 {
41 unsigned int old_cw;
42 fpu_fix_start(&old_cw);
43
44 typedef GeMatrix<FullStorage<T> > Matrix;
45 typedef DenseVector<Array<T> > Vector;
46 typedef Matrix::IndexType IndexType;
47 typedef DenseVector<Array<IndexType> > IndexVector;
48
49 const Underscore<IndexType> _;
50
51 const IndexType m = 4,
52 n = 5;
53
54 Matrix Ab(m, n);
55 IndexVector piv(m);
56
57 Ab = 2, 3, -1, 0, 20,
58 -6, -5, 0, 2, -33,
59 2, -5, 6, -6, -43,
60 4, 6, 2, -3, 49;
61
62 cout << "Ab = " << Ab << endl;
63
64 lapack::trf(Ab, piv);
65
66 const auto A = Ab(_,_(1,m));
67 auto B = Ab(_,_(m+1,n));
68
69 blas::sm(Left, NoTrans, T(1), A.upper(), B);
70
71 cout << "X = " << B << endl;
72
73 ///
74 /// ... and at the end restore FPU rounding behavior as mentioned above.
75 ///
76 fpu_fix_end(&old_cw);
77 }
78
79 ///
80 /// :links: __QD Library__ -> http://crd-legacy.lbl.gov/~dhbailey/mpdist/
81 ///