Macaulay2 Engine
Loading...
Searching...
No Matches
dmat-ffpack.cpp
Go to the documentation of this file.
1// Copyright 2005-2012 Michael E. Stillman
2
30
31#if 0
32// This file is not in use. These functions are now in dmat.cpp.
33#if 0
34 template<typename CoeffRing >
35 template<class RingType>
36 size_t DMat < CoeffRing >::rank(typename enable_if<is_givaro_or_ffpack<RingType>::value >::type* dummy ) const
37 {
38 // assert not necessary because the test is already done by "enable_if<is_givaro_or_ffpack<RingType>::value >"
39 // assert( typeid(CoeffRing) == typeid(M2::ARingZZpFFPACK) || typeid(CoeffRing) == typeid(M2::ARingGFGivaro ));
40 std::cout << "Calling rankGF_or_FFPACK" << std::endl;
41 ElementType *N = newarray(ElementType, n_rows() * n_cols() );
45 copy_elems( n_rows()*n_cols(), N, 1, get_array(), 1);
48
49 /* //debug
50 typename MatrixType::ElementType *Npos=N;
51 for ( int currRow=0; currRow < n_rows(); currRow++ )
52 for ( int currCol =0; currCol < n_cols(); currCol++ )
53 {
54 typename MatrixType::ElementType entry;
55 get_entry(currRow, currCol,entry) ;
56 ring().field().init(Npos, entry );
57 // mat.setEntry( currRow, currCol ,( (int)rand() ) % characteristic );
58 }*/
59
60 size_t result = FFPACK::Rank(ring().field(), n_cols(), n_rows(), N, n_rows() );
61 freemem(N);
62 return result;
63 }
64#endif
65
66 template<typename CoeffRing>
67 size_t FFpackRank(const DMat<CoeffRing>& mat)
68 {
69 typedef typename CoeffRing::ElementType ElementType;
70 std::cout << "Calling FFpackRank" << std::endl;
71 assert( typeid(CoeffRing) == typeid(M2::ARingZZpFFPACK) || typeid(CoeffRing) == typeid(M2::ARingGFGivaro ));
72 ElementType* N = newarray( ElementType, mat.n_rows() * mat.n_cols());
73 mat.copy_elems(mat.n_rows()*mat.n_cols(), N, 1, mat.get_array(), 1);
76 size_t result = FFPACK::Rank(mat.ring().field(), mat.n_cols(), mat.n_rows(), N, mat.n_rows());
77 freemem(N);
78 return result;
79 }
80
81 template<typename CoeffRing>
82 void FFpackDeterminant(const DMat<CoeffRing>& mat,
83 typename CoeffRing::ElementType& result )
84 {
85 typedef typename CoeffRing::ElementType ElementType;
86 std::cout << "Calling FFpackDeterminant" << std::endl;
87 assert( typeid(CoeffRing) == typeid(M2::ARingZZpFFPACK) || typeid(CoeffRing) == typeid(M2::ARingGFGivaro ));
88 ElementType* N = newarray( ElementType, mat.n_rows() * mat.n_cols());
89 mat.copy_elems(mat.n_rows()*mat.n_cols(), N, 1, mat.get_array(), 1);
92 result = FFPACK::Det(mat.ring().field(), mat.n_cols(), mat.n_rows(), N, mat.n_rows());
93 freemem(N);
94 }
95
96 template<typename CoeffRing>
97 bool FFpackInvert(const DMat<CoeffRing> &mat, DMat<CoeffRing> &inverse)
98 {
99 typedef typename CoeffRing::ElementType ElementType;
100 assert(mat.n_rows() == mat.n_cols());
101 ElementType* N = newarray( ElementType, mat.n_rows() * mat.n_cols());
102 mat.copy_elems(mat.n_rows()*mat.n_cols(), N, 1, mat.get_array(), 1);
103
104 size_t n = mat.n_rows(); // same as n_cols()
105 int nullspacedim;
106 FFPACK::Invert2(mat.ring().field(), n, N, n, inverse.get_array(), n, nullspacedim);
107
108 freemem(N);
109 return true;
110 }
111
112 template<typename CoeffRing>
113 void FFpackNullSpace(const DMat<CoeffRing> &mat,
114 DMat<CoeffRing> &nullspace,
115 bool right_side)
116 {
117 right_side = !right_side; // because FFPACK stores by rows, not by columns.
118
119 typedef typename CoeffRing::ElementType ElementType;
120 ElementType* N = newarray( ElementType, mat.n_rows() * mat.n_cols());
121 mat.copy_elems(mat.n_rows()*mat.n_cols(), N, 1, mat.get_array(), 1);
122
123 size_t nr = mat.n_rows();
124 size_t nc = mat.n_cols();
125
126 ElementType *nullspaceFFPACK = 0;
127
128 size_t nullspace_dim;
129 size_t nullspace_leading_dim;
130
131 FFPACK::NullSpaceBasis(mat.ring().field(),
132 (right_side ? FFLAS::FflasRight : FFLAS::FflasLeft),
133 nc, nr, N, nr, nullspaceFFPACK, nullspace_leading_dim, nullspace_dim);
134
135 std::cerr << "leading dim = " << nullspace_leading_dim << " and dim = " << nullspace_dim << std::endl;
136 //NOTUSED? size_t nullspace_nrows = (right_side ? nc : nullspace_dim);
137 if (right_side && nullspace_dim != nullspace_leading_dim)
138 {
139 std::cerr << "error: this should not happen!" << std::endl;
140 }
141 else if (!right_side && nullspace_leading_dim != nc)
142 {
143 std::cerr << "error: this should not happen either!" << std::endl;
144 }
145
146 if (right_side)
147 nullspace.resize(nullspace_dim,nr);
148 else
149 nullspace.resize(nc,nullspace_dim);
150
151 mat.copy_elems(nullspace.n_rows() * nullspace.n_cols(), nullspace.get_array(), 1, nullspaceFFPACK, 1);
152
153 delete [] nullspaceFFPACK;
154 }
155
156 template<typename CoeffRing>
157 M2_arrayintOrNull FFpackRankProfile(const DMat<CoeffRing> &mat,
158 bool row_profile)
159 {
160 // Note that FFPack stores matrices by row, not column, the opposite of what we do.
161 // So row_profile true means use ffpack column rank profile!
162
163 typedef typename CoeffRing::ElementType ElementType;
164 ElementType* N = newarray( ElementType, mat.n_rows() * mat.n_cols());
165 mat.copy_elems(mat.n_rows()*mat.n_cols(), N, 1, mat.get_array(), 1);
166
167 size_t * prof; // this is where the result will be placed
168
169 size_t rk;
170 if (!row_profile)
171 rk = FFPACK::RowRankProfile(mat.ring().field(),
172 mat.n_cols(),mat.n_rows(),
173 N,mat.n_rows(),
174 prof);
175 else
176 rk = FFPACK::ColumnRankProfile(mat.ring().field(),
177 mat.n_cols(),mat.n_rows(),
178 N,mat.n_rows(),
179 prof);
180
181 M2_arrayint profile = M2_makearrayint(static_cast<int>(rk));
182 for (size_t i=0; i<rk; i++)
183 profile->array[i] = static_cast<int>(prof[i]);
184
185 delete [] prof;
186 freemem(N);
187
188 return profile;
189 }
190
191 template<typename CoeffRing>
192 bool FFpackSolveLinear(const DMat<CoeffRing> &mat,
193 DMat<CoeffRing> &X,
194 const DMat<CoeffRing> &B,
195 bool right_side)
196 {
197 std::cerr << "inside FFpackSolveLinear" << std::endl;
198
199 typedef typename CoeffRing::ElementType ElementType;
200 size_t a_rows = mat.n_rows();
201 size_t a_cols = mat.n_cols();
202
203 size_t b_rows = B.n_rows();
204 size_t b_cols = B.n_cols();
205
206 ElementType* ffpackA = newarray(ElementType, mat.n_rows() * mat.n_cols());
207 mat.copy_elems(mat.n_rows()*mat.n_cols(), ffpackA, 1, mat.get_array(), 1);
208
209 ElementType* ffpackB = newarray(ElementType, b_rows * b_cols);
210 B.copy_elems(b_rows * b_cols, ffpackB, 1, B.get_array(), 1);
211
212 // preallocate the space for the solutions:
213 size_t x_rows = (right_side ? a_cols : b_rows);
214 size_t x_cols = (right_side ? b_cols : a_rows);
215 //NOTUSED? size_t n_eqns = (right_side ? b_cols : b_rows);
216
217 ElementType *ffpackX = newarray_clear(ElementType, x_rows * x_cols);
218
219 int info; // >0 if the system is inconsistent, ==0 means success
220
221 FFPACK::fgesv(mat.ring().field(),
222 (!right_side ? FFLAS::FflasLeft : FFLAS::FflasRight),
223 a_cols, a_rows,
224 (!right_side ? b_cols : b_rows),
225 ffpackA,
226 a_rows, // leading dim of A
227 ffpackX, x_rows,
228 ffpackB, b_rows,
229 &info);
230
231 if (info > 0)
232 {
233 // the system is inconsistent
234 ERROR("the system is inconsistent");
235 return false;
236 }
237
238 X.resize(x_rows, x_cols);
239 X.copy_elems(x_rows * x_cols, X.get_array(), 1, ffpackX, 1);
240
241 delete [] ffpackX;
242
243 return true;
244 }
245
246
247
248
249 template<typename CoeffRing>
250 void FFpackScalarMultiplyInPlace(DMat<CoeffRing>& A,
251
252 const typename CoeffRing::ElementType& scalar)
253 {
254 size_t incx = 1;
255 FFLAS::fgemm( A.ring().field(), A.n_rows()*A.n_cols(), scalar, A.get_array(), incx);
256 }
257
258
259 template<typename CoeffRing>
260 void FFpackMatrixAdd(DMat<CoeffRing>& C,
261 const DMat<CoeffRing>& A,
262 const DMat<CoeffRing>& B)
263 {
264
265 FFLAS::fadd ( C.ring().field(),
266 C.n_cols(), C.n_rows(),
267 A.get_array(), A.n_rows(),
268 B.get_array(), B.n_rows(),
269 C.get_array(), C.n_rows() );
270 }
271
272 template<typename CoeffRing>
273 void FFpackMatrixAddInPlace(DMat<CoeffRing>& C,
274 const DMat<CoeffRing>& A)
275 {
276
277 FFLAS::faddin ( C.ring().field(),
278 C.n_cols(), C.n_rows(),
279 A.get_array(), A.n_rows(),
280 C.get_array(), C.n_rows() );
281 }
282
283 template<typename CoeffRing>
284 void FFpackMatrixSub(DMat<CoeffRing>& C,
285 const DMat<CoeffRing>& A,
286 const DMat<CoeffRing>& B)
287 {
288
289 FFLAS::fsub ( C.ring().field(),
290 C.n_cols(), C.n_rows(),
291 A.get_array(), A.n_rows(),
292 B.get_array(), B.n_rows(),
293 C.get_array(), C.n_rows() );
294 }
295
296 // there is no BLAS or even FFPACK routine for transposing a matrix.
297
298 template<typename CoeffRing>
299 void FFpackAddMultipleTo(DMat<CoeffRing>& C,
300 const DMat<CoeffRing>& A,
301 const DMat<CoeffRing>& B,
302 bool transposeA,
303 bool transposeB,
304 const typename CoeffRing::ElementType& a,
305 const typename CoeffRing::ElementType& b)
306 /* A,B,C should be mutable matrices over a finite prime field, and a,b
307 elements of this field.
308 C = b*C + a * op(A)*op(B),
309 where op(A) = A or transpose(A), depending on transposeA
310 where op(B) = B or transpose(B), depending on transposeB
311 */
312 {
313
314 FFLAS::FFLAS_TRANSPOSE tA = (transposeA ? FFLAS::FflasTrans : FFLAS::FflasNoTrans);
315 FFLAS::FFLAS_TRANSPOSE tB = (transposeB ? FFLAS::FflasTrans : FFLAS::FflasNoTrans);
316
317 size_t m = (transposeB ? B.n_rows() : B.n_cols());
318 size_t n = (transposeA ? A.n_cols() : A.n_rows());
319
320 size_t k = (transposeA ? A.n_rows() : A.n_cols());
321 size_t k2 = (transposeB ? B.n_cols() : B.n_rows());
322
323 assert(k == k2); // The user of this function must insure that sizes are correct.
324 if (k!=k2)
325 {
326 ERROR("matrices are not composable. this error is handled in Macaulay2 code and this message should never appear! ");
327 return;
328 }
329
330
331 FFLAS::fgemm( C.ring().field(),
332 tB, tA,
333 m,n,k,
334 a,
335 B.get_array(),
336 B.n_rows(),
337 A.get_array(),
338 A.n_rows(),
339 b,
340 C.get_array(),
341 C.n_rows()
342 );
343 }
344
346 // ARingZZpFFPACK specific linear algebra functions //
348 template<>
349 void DMat<M2::ARingZZpFFPACK>::addInPlace(const DMat& B)
350 {
351 std::cout << "DMat<M2::ARingZZpFFPACK>::addInPlace()" << std::endl;
352 FFpackMatrixAddInPlace<M2::ARingZZpFFPACK>(*this,B);
353 }
354
355
356 template<>
357 size_t DMat<M2::ARingZZpFFPACK>::rank() const
358 {
359 std::cout << "DMat<M2::ARingZZpFFPACK>::rank()" << std::endl;
360 return FFpackRank<M2::ARingZZpFFPACK>(*this);
361 }
362
363 template<>
364 void DMat<M2::ARingZZpFFPACK>::determinant(elem &result) const
365 {
366 std::cout << "Calling DMat<M2::ARingZZpFFPACK>::determinant" << std::endl;
367 FFpackDeterminant<M2::ARingZZpFFPACK>(*this, result);
368 }
369
370 template<>
371 bool DMat<M2::ARingZZpFFPACK>::invert(DMat<M2::ARingZZpFFPACK> &inverse) const
372 {
373 std::cout << "Calling DMat<M2::ARingZZpFFPACK>::inverse" << std::endl;
374 return FFpackInvert<M2::ARingZZpFFPACK>(*this, inverse);
375 }
376
377 template<>
378 void DMat<M2::ARingZZpFFPACK>::nullSpace(DMat<M2::ARingZZpFFPACK> &nullspace, bool right_side) const
379 {
380 std::cout << "Calling DMat<M2::ARingZZpFFPACK>::nullspace" << std::endl;
381 FFpackNullSpace<M2::ARingZZpFFPACK>(*this, nullspace, right_side);
382 }
383
384 template<>
385 M2_arrayintOrNull DMat<M2::ARingZZpFFPACK>::rankProfile(bool row_profile) const
386 {
387 std::cout << "Calling DMat<M2::ARingZZpFFPACK>::rankProfile" << std::endl;
388 return FFpackRankProfile(*this, row_profile);
389 }
390
391 template<>
392 bool DMat<M2::ARingZZpFFPACK>::solveLinear(DMat<M2::ARingZZpFFPACK> &X,
393 const DMat<M2::ARingZZpFFPACK> &B,
394 bool right_side) const
395 {
396 std::cout << "Calling DMat<M2::ARingZZpFFPACK>::solveLinear" << std::endl;
397 return FFpackSolveLinear(*this, X, B, right_side);
398 }
399
400 template<>
401 void DMat<M2::ARingZZpFFPACK>::addMultipleTo(const DMat<M2::ARingZZpFFPACK> &A,
402 const DMat<M2::ARingZZpFFPACK> &B,
403 bool transposeA,
404 bool transposeB,
405 const ElementType& a,
406 const ElementType& b)
407 {
408 std::cout << "Calling DMat<M2::ARingZZpFFPACK>::addMultipleTo *" << std::endl;
409 FFpackAddMultipleTo(*this, A, B, transposeA, transposeB, a, b);
410 }
411
412#if 1
413
415 // ARingGFGivaro specific linear algebra functions /////////
417
418 template<>
419 void DMat<M2::ARingGFGivaro>::addInPlace(const DMat& B)
420 {
421 std::cout << "not implemented yet" << std::endl;
422 }
423
424 template<>
425 size_t DMat<M2::ARingGFGivaro>::rank() const
426 {
427 std::cout << "Calling DMat<M2::ARingGFGivaro>::rank()" << std::endl;
428 return FFpackRank<M2::ARingGFGivaro>(*this);
429 }
430
431 template<>
432 void DMat<M2::ARingGFGivaro>::determinant(elem &result) const
433 {
434 std::cout << "Calling DMat<M2::ARingGFGivaro>::determinant" << std::endl;
435 FFpackDeterminant<M2::ARingGFGivaro>(*this, result );
436 }
437
438 template<>
439 bool DMat<M2::ARingGFGivaro>::invert(DMat<M2::ARingGFGivaro> &inverse) const
440 {
441 std::cout << "Calling DMat<M2::ARingGFGivaro>::inverse" << std::endl;
442 return FFpackInvert<M2::ARingGFGivaro>(*this, inverse);
443 }
444
445 template<>
446 void DMat<M2::ARingGFGivaro>::nullSpace(DMat<M2::ARingGFGivaro> &nullspace, bool right_side) const
447 {
448 std::cout << "Calling DMat<M2::ARingGFGivaro>::nullspace" << std::endl;
449 FFpackNullSpace<M2::ARingGFGivaro>(*this, nullspace, right_side);
450 }
451
452 template<>
453 M2_arrayintOrNull DMat<M2::ARingGFGivaro>::rankProfile(bool row_profile) const
454 {
455 std::cout << "Calling DMat<M2::ARingGFGivaro>::rankProfile" << std::endl;
456 return FFpackRankProfile(*this, row_profile);
457 }
458
459 template<>
460 bool DMat<M2::ARingGFGivaro>::solveLinear(DMat<M2::ARingGFGivaro> &X,
461 const DMat<M2::ARingGFGivaro> &B,
462 bool right_side) const
463 {
464 std::cout << "Calling DMat<M2::ARingGFGivaro>::solveLinear" << std::endl;
465 return FFpackSolveLinear(*this, X, B, right_side);
466 }
467
468 template<>
469 void DMat<M2::ARingGFGivaro>::addMultipleTo(const DMat<M2::ARingGFGivaro> &A,
470 const DMat<M2::ARingGFGivaro> &B,
471 bool transposeA,
472 bool transposeB,
473 const ElementType& a,
474 const ElementType& b)
475 {
476 std::cout << "Calling DMat<M2::ARingGFGivaro>::addMultipleTo" << std::endl;
477 FFpackAddMultipleTo(*this, A, B, transposeA, transposeB, a, b);
478 }
479
480#endif
481#endif
482
483// Local Variables:
484// compile-command: "make -C $M2BUILDDIR/Macaulay2/e "
485// indent-tabs-mode: nil
486// End:
DMat::resize
void resize(size_t new_nrows, size_t new_ncols)
Definition dmat.hpp:157
DMat::ring
const ACoeffRing & ring() const
Definition dmat.hpp:143
DMat
Definition dmat.hpp:62
M2::ARingZZpFFPACK
wrapper for the FFPACK::ModularBalanced<double> field implementation
Definition aring-zzp-ffpack.hpp:64
M2_arrayint
freemem
void freemem(void *s)
Definition m2-mem.cpp:103
ERROR
const int ERROR
Definition m2-mem.cpp:55
result
VALGRIND_MAKE_MEM_DEFINED & result(result)
M2_makearrayint
M2_arrayint M2_makearrayint(int n)
Definition m2-types.cpp:6
M2_arrayintOrNull
M2_arrayint M2_arrayintOrNull
Definition m2-types.h:99
ffpackInterface::inverse
bool inverse(const DMatZZpFFPACK &A, DMatZZpFFPACK &result_inv)
Definition dmat.cpp:113
newarray
#define newarray(T, len)
Definition newdelete.hpp:82
newarray_clear
#define newarray_clear(T, len)
Definition newdelete.hpp:83