1 #include <iostream>
2
3
4 /*
5 #define CXXBLAS_DEBUG_OUT(X) \
6 flens::verbose::ClosureLog::append(false) << " => BLAS: " << X << ";"; \
7 */
8
9
10 #include <flens/debug/aux/aux.h>
11 #include <flens/flens.cxx>
12 #include <flens/debug/aux/aux.tcc>
13
14 using namespace std;
15 using namespace flens;
16
17 typedef double T;
18
19 int
20 main()
21 {
22 typedef GeMatrix<FullStorage<T> > GeMatrix;
23 typedef DenseVector<Array<T> > DenseVector;
24 typedef DenseVector::IndexType IndexType;
25
26 const Underscore<IndexType> _;
27
28 int n = 3;
29
30 verbose::ClosureLog::start("log");
31
32 GeMatrix A_(n,n), B_(n,n), C_(n,n), D_(n,n), E_(n,n);
33
34 GeMatrix::View A = A_(_,_(1,3)),
35 B = B_(_,_(1,3)),
36 C = C_(_,_(1,3)),
37 D = D_(_,_(1,3)),
38 E = E_(_,_(1,3));
39
40 A = 1;
41 B = 2;
42 C = 3;
43 D = 4;
44 E = 5;
45
46 DenseVector a(n), b(n), c(n), d(n), e(n);
47
48 a = 1;
49 b = 2;
50 c = 3;
51 d = 4;
52 e = 5;
53
54 /*
55 b = A.upper()*b;
56 c = A.upper()*b;
57 c += A.upper()*b;
58
59 C = transpose(A*B);
60 C = A*B;
61 A = A*A;
62 C = C + transpose(A)*B;
63 C = 2*C + transpose(A)*B;
64 C = 2*C + transpose(transpose(A)*B);
65 C = transpose(2*C + transpose(A)*B);
66
67 b = 2*b + 3*transpose(A)*a;
68 b = 3*transpose(A)*a + 2*b;
69 */
70
71 cout << "HasEngine<ScalarValue<int> >::value = "
72 << HasEngine<ScalarValue<int> >::value
73 << endl;
74
75 cout << "HasEngine<GeMatrix>::value = "
76 << HasEngine<GeMatrix>::value
77 << endl;
78
79 cout << "HasEngine<GeMatrix::TrMatrixView>::value = "
80 << HasEngine<GeMatrix::TriangularView>::value
81 << endl;
82
83 cout << "IsFullStorage<GeMatrix::Engine>::value = "
84 << IsFullStorage<GeMatrix::Engine>::value
85 << endl;
86
87 cout << "HasFullStorage<GeMatrix>::value = "
88 << HasFullStorage<GeMatrix>::value
89 << endl;
90
91 cout << "IsFullStorage<int>::value = "
92 << IsFullStorage<int>::value
93 << endl;
94
95 using DEBUGCLOSURE::identical;
96 cout << "identical(C, C.upper()) = " << identical(C, C.upper()) << endl;
97
98 C = A*C.upper();
99 // C = C.upper()*A;
100 // C = A.upper()*C.upper();
101
102 /*
103 C = A.upper();
104
105 a = c - d;
106
107 c = b - A*a;
108 c = b - 2*A*a;
109 c = 2*c + transpose(2*transpose(A))*a;
110 c = 2*c + transpose(2*transpose(A+B))*(a+b);
111
112 e = 2*b + A*(a+c);
113 cout << "e = " << e << endl;
114
115 e = 2*b + transpose(A)*(a+c);
116 cout << "e = " << e << endl;
117
118 e = 2*b + transpose(3*A)*(a+c);
119 cout << "e = " << e << endl;
120
121 e = 2*b + transpose(3*A + transpose(A))*(a+c);
122 cout << "e = " << e << endl;
123
124
125 cout << "A = " << A << endl;
126 A += transpose(A);
127 cout << "A = " << A << endl;
128
129 A += transpose(A);
130 cout << "A = " << A << endl;
131
132 GeMatrix X(n, n);
133 X = A;
134 cout << "X = " << X << endl;
135
136 X(_(1,2),_).upper() = transpose(E(_,_(1,2)).lower());
137 cout << "X = " << X << endl;
138
139
140 b = A*b;
141
142 b = A.upper()*b;
143 c += 3*(2*transpose(A.upper())*b);
144 b += A.upper()*b;
145
146 B = C.lower();
147 cout << "A = " << A << endl;
148 cout << "C = " << C << endl;
149 cout << "A.upper() = " << A.upper() << endl;
150 cout << "C.lower() = " << C.lower() << endl;
151 cout << "B = " << B << endl;
152
153 D = C.lower();
154 cout << "D = " << D << endl;
155
156
157 A = 3*B + (4*C + 5*D);
158 cerr << "A = " << A << endl;
159
160 A = 3*B + (4*C + 2*(5*D));
161 cerr << "A = " << A << endl;
162
163 A = 3*B + (4*C + D/5);
164 cerr << "A = " << A << endl;
165
166 A = 3*B + (4*C + D/5/2);
167 cerr << "A = " << A << endl;
168
169 A = B + (C + D);
170 cerr << "A = " << A << endl;
171
172 A = B + (C - D);
173 cerr << "A = " << A << endl;
174
175 B(1,_) = 4;
176 cerr << "B = " << B << endl;
177
178 A = transpose(B);
179 cerr << "A = " << A << endl;
180
181 A = transpose(A);
182 cerr << "A = " << A << endl;
183
184 A = transpose(A(_(1,n),_(1,n)));
185 cerr << "A = " << A << endl;
186
187 A = 3*B;
188 cerr << "A = " << A << endl;
189
190 A = 3*A;
191 cerr << "A = " << A << endl;
192
193 A = 2*(3*B);
194 cerr << "A = " << A << endl;
195
196 A = B/3;
197 cerr << "A = " << A << endl;
198
199 A = A/3;
200 cerr << "A = " << A << endl;
201
202 A = B/3/2;
203 cerr << "A = " << A << endl;
204
205
206 a = a + b;
207
208 a = b + a;
209
210 a = a + c + d + e+ b + a(_(1,n));
211
212 a = a + a(_(1,n)) + d + e+ b + a(_(1,n));
213
214 a = 3*a;
215
216 a = a*3*2;
217
218 e = 2*a - b/2 + 2*c + 5*d;
219 cerr << "e = " << e << endl;
220
221 e = M_PI*(2*a) - b/2 + 2*(3*c) + d;
222 cerr << "e = " << e << endl;
223
224 e = M_PI*(2*a) - b/2/3 + 2*(3*c) + d;
225 cerr << "e = " << e << endl;
226
227 e = b/2;
228 cerr << "e = " << e << endl;
229
230 const DenseVector &x = a;
231 cerr << "x = " << x << endl;
232
233 const DenseVector &y = b;
234 cerr << "y = " << y << endl;
235
236 e = 2*x*y*x;
237 cerr << "e = " << e << endl;
238
239 e += (x+x)*(y*x);
240 cerr << "e = " << e << endl;
241
242 e = ((a - b) - c) - d;
243 cerr << "e = " << e << endl;
244
245 e = a - (b - (c + d));
246 cerr << "e = " << e << endl;
247
248 e = a + (b + (c + d));
249 cerr << "e = " << e << endl;
250
251
252 D = transpose(C.lower())*transpose(A+B);
253
254 cout << "D = " << D << endl;
255 cout << "C = " << C << endl;
256 cout << "C.lower().symmetric() = " << C.lower().symmetric() << endl;
257 GeMatrix XX = C.lower().symmetric() + C;
258 cout << "-> XX = " << XX << endl;
259
260 D = C*A;
261 cout << "D = " << D << endl;
262 D = C.lower().symmetric()*D;
263 cout << "D = " << D << endl;
264 D = D.lower().symmetric()*A;
265
266 cout << "D = " << D << endl;
267 b = (D.lower().symmetric()+D)*a;
268 b = D.lower().symmetric()*a;
269
270 A = E.lower() + C.upper();
271 cout << "E.lower() = " << E.lower() << endl;
272 cout << "C.upper() = " << C.upper() << endl;
273 cout << "A = " << A << endl;
274 */
275
276 verbose::ClosureLog::stop();
277 }
2
3
4 /*
5 #define CXXBLAS_DEBUG_OUT(X) \
6 flens::verbose::ClosureLog::append(false) << " => BLAS: " << X << ";"; \
7 */
8
9
10 #include <flens/debug/aux/aux.h>
11 #include <flens/flens.cxx>
12 #include <flens/debug/aux/aux.tcc>
13
14 using namespace std;
15 using namespace flens;
16
17 typedef double T;
18
19 int
20 main()
21 {
22 typedef GeMatrix<FullStorage<T> > GeMatrix;
23 typedef DenseVector<Array<T> > DenseVector;
24 typedef DenseVector::IndexType IndexType;
25
26 const Underscore<IndexType> _;
27
28 int n = 3;
29
30 verbose::ClosureLog::start("log");
31
32 GeMatrix A_(n,n), B_(n,n), C_(n,n), D_(n,n), E_(n,n);
33
34 GeMatrix::View A = A_(_,_(1,3)),
35 B = B_(_,_(1,3)),
36 C = C_(_,_(1,3)),
37 D = D_(_,_(1,3)),
38 E = E_(_,_(1,3));
39
40 A = 1;
41 B = 2;
42 C = 3;
43 D = 4;
44 E = 5;
45
46 DenseVector a(n), b(n), c(n), d(n), e(n);
47
48 a = 1;
49 b = 2;
50 c = 3;
51 d = 4;
52 e = 5;
53
54 /*
55 b = A.upper()*b;
56 c = A.upper()*b;
57 c += A.upper()*b;
58
59 C = transpose(A*B);
60 C = A*B;
61 A = A*A;
62 C = C + transpose(A)*B;
63 C = 2*C + transpose(A)*B;
64 C = 2*C + transpose(transpose(A)*B);
65 C = transpose(2*C + transpose(A)*B);
66
67 b = 2*b + 3*transpose(A)*a;
68 b = 3*transpose(A)*a + 2*b;
69 */
70
71 cout << "HasEngine<ScalarValue<int> >::value = "
72 << HasEngine<ScalarValue<int> >::value
73 << endl;
74
75 cout << "HasEngine<GeMatrix>::value = "
76 << HasEngine<GeMatrix>::value
77 << endl;
78
79 cout << "HasEngine<GeMatrix::TrMatrixView>::value = "
80 << HasEngine<GeMatrix::TriangularView>::value
81 << endl;
82
83 cout << "IsFullStorage<GeMatrix::Engine>::value = "
84 << IsFullStorage<GeMatrix::Engine>::value
85 << endl;
86
87 cout << "HasFullStorage<GeMatrix>::value = "
88 << HasFullStorage<GeMatrix>::value
89 << endl;
90
91 cout << "IsFullStorage<int>::value = "
92 << IsFullStorage<int>::value
93 << endl;
94
95 using DEBUGCLOSURE::identical;
96 cout << "identical(C, C.upper()) = " << identical(C, C.upper()) << endl;
97
98 C = A*C.upper();
99 // C = C.upper()*A;
100 // C = A.upper()*C.upper();
101
102 /*
103 C = A.upper();
104
105 a = c - d;
106
107 c = b - A*a;
108 c = b - 2*A*a;
109 c = 2*c + transpose(2*transpose(A))*a;
110 c = 2*c + transpose(2*transpose(A+B))*(a+b);
111
112 e = 2*b + A*(a+c);
113 cout << "e = " << e << endl;
114
115 e = 2*b + transpose(A)*(a+c);
116 cout << "e = " << e << endl;
117
118 e = 2*b + transpose(3*A)*(a+c);
119 cout << "e = " << e << endl;
120
121 e = 2*b + transpose(3*A + transpose(A))*(a+c);
122 cout << "e = " << e << endl;
123
124
125 cout << "A = " << A << endl;
126 A += transpose(A);
127 cout << "A = " << A << endl;
128
129 A += transpose(A);
130 cout << "A = " << A << endl;
131
132 GeMatrix X(n, n);
133 X = A;
134 cout << "X = " << X << endl;
135
136 X(_(1,2),_).upper() = transpose(E(_,_(1,2)).lower());
137 cout << "X = " << X << endl;
138
139
140 b = A*b;
141
142 b = A.upper()*b;
143 c += 3*(2*transpose(A.upper())*b);
144 b += A.upper()*b;
145
146 B = C.lower();
147 cout << "A = " << A << endl;
148 cout << "C = " << C << endl;
149 cout << "A.upper() = " << A.upper() << endl;
150 cout << "C.lower() = " << C.lower() << endl;
151 cout << "B = " << B << endl;
152
153 D = C.lower();
154 cout << "D = " << D << endl;
155
156
157 A = 3*B + (4*C + 5*D);
158 cerr << "A = " << A << endl;
159
160 A = 3*B + (4*C + 2*(5*D));
161 cerr << "A = " << A << endl;
162
163 A = 3*B + (4*C + D/5);
164 cerr << "A = " << A << endl;
165
166 A = 3*B + (4*C + D/5/2);
167 cerr << "A = " << A << endl;
168
169 A = B + (C + D);
170 cerr << "A = " << A << endl;
171
172 A = B + (C - D);
173 cerr << "A = " << A << endl;
174
175 B(1,_) = 4;
176 cerr << "B = " << B << endl;
177
178 A = transpose(B);
179 cerr << "A = " << A << endl;
180
181 A = transpose(A);
182 cerr << "A = " << A << endl;
183
184 A = transpose(A(_(1,n),_(1,n)));
185 cerr << "A = " << A << endl;
186
187 A = 3*B;
188 cerr << "A = " << A << endl;
189
190 A = 3*A;
191 cerr << "A = " << A << endl;
192
193 A = 2*(3*B);
194 cerr << "A = " << A << endl;
195
196 A = B/3;
197 cerr << "A = " << A << endl;
198
199 A = A/3;
200 cerr << "A = " << A << endl;
201
202 A = B/3/2;
203 cerr << "A = " << A << endl;
204
205
206 a = a + b;
207
208 a = b + a;
209
210 a = a + c + d + e+ b + a(_(1,n));
211
212 a = a + a(_(1,n)) + d + e+ b + a(_(1,n));
213
214 a = 3*a;
215
216 a = a*3*2;
217
218 e = 2*a - b/2 + 2*c + 5*d;
219 cerr << "e = " << e << endl;
220
221 e = M_PI*(2*a) - b/2 + 2*(3*c) + d;
222 cerr << "e = " << e << endl;
223
224 e = M_PI*(2*a) - b/2/3 + 2*(3*c) + d;
225 cerr << "e = " << e << endl;
226
227 e = b/2;
228 cerr << "e = " << e << endl;
229
230 const DenseVector &x = a;
231 cerr << "x = " << x << endl;
232
233 const DenseVector &y = b;
234 cerr << "y = " << y << endl;
235
236 e = 2*x*y*x;
237 cerr << "e = " << e << endl;
238
239 e += (x+x)*(y*x);
240 cerr << "e = " << e << endl;
241
242 e = ((a - b) - c) - d;
243 cerr << "e = " << e << endl;
244
245 e = a - (b - (c + d));
246 cerr << "e = " << e << endl;
247
248 e = a + (b + (c + d));
249 cerr << "e = " << e << endl;
250
251
252 D = transpose(C.lower())*transpose(A+B);
253
254 cout << "D = " << D << endl;
255 cout << "C = " << C << endl;
256 cout << "C.lower().symmetric() = " << C.lower().symmetric() << endl;
257 GeMatrix XX = C.lower().symmetric() + C;
258 cout << "-> XX = " << XX << endl;
259
260 D = C*A;
261 cout << "D = " << D << endl;
262 D = C.lower().symmetric()*D;
263 cout << "D = " << D << endl;
264 D = D.lower().symmetric()*A;
265
266 cout << "D = " << D << endl;
267 b = (D.lower().symmetric()+D)*a;
268 b = D.lower().symmetric()*a;
269
270 A = E.lower() + C.upper();
271 cout << "E.lower() = " << E.lower() << endl;
272 cout << "C.upper() = " << C.upper() << endl;
273 cout << "A = " << A << endl;
274 */
275
276 verbose::ClosureLog::stop();
277 }