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, ColMajor> >   Matrix;
45     typedef DenseVector<Array<T> >                Vector;
46 
47     const int n = 5;
48 
49     Matrix   A(n, n), VL(n, n), VR(n, n);
50     Vector   wr(n), wi(n);
51 
52 
53     A =  2,   3,  -1,   0,  2,
54         -6,  -5,   0,   2, -6,
55          2,  -5,   6,  -6,  2,
56          2,   3,  -1,   0,  8,
57         -6,  -5,  10,   2, -6;
58 
59     cerr << "A = " << A << endl;
60 
61     Vector   work;
62     lapack::ev(truetrue, A, wr, wi, VL, VR, work);
63 
64     cerr << "wr = " << wr << endl;
65     cerr << "wi = " << wi << endl;
66     cerr << "VL = " << VL << endl;
67     cerr << "VR = " << VR << endl;
68 
69     ///
70     ///  ... and at the end restore FPU rounding behavior as mentioned above.
71     ///
72     fpu_fix_end(&old_cw);
73 }
74 
75 ///
76 ///  :links: __QD Library__ -> http://crd-legacy.lbl.gov/~dhbailey/mpdist/
77 ///