1 /* 2 * Copyright (C) 2015 FFLAS-FFPACK 3 * 4 * Written by Brice Boyer (briceboyer) <boyer.brice@gmail.com> 5 * 6 * 7 * ========LICENCE======== 8 * This file is part of the library FFLAS-FFPACK. 9 * 10 * FFLAS-FFPACK is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 23 * ========LICENCE======== 24 *. 25 */ 26 27 /** @file ffpack-c.h 28 * @author Brice Boyer 29 * @brief C functions calls for FFPACK 30 * @see ffpack/ffpack.h 31 */ 32 33 #ifndef __FFLASFFPACK_interfaces_libs_ffpack_c_H 34 #define __FFLASFFPACK_interfaces_libs_ffpack_c_H 35 //#include "fflas-ffpack/fflas-ffpack-config.h" 36 37 #ifndef FFPACK_COMPILED 38 #define FFPACK_COMPILED 39 #endif 40 41 #include <stdbool.h> 42 #include <stdlib.h> 43 #include <inttypes.h> 44 45 #ifdef __cplusplus 46 extern "C" { 47 #endif 48 49 50 #ifndef __FFLASFFPACK_interfaces_libs_fflas_c_H 51 52 enum FFLAS_C_ORDER { 53 FflasRowMajor=101, 54 FflasColMajor=102 55 }; 56 enum FFLAS_C_TRANSPOSE { 57 FflasNoTrans = 111, 58 FflasTrans = 112 59 }; 60 enum FFLAS_C_UPLO { 61 FflasUpper = 121, 62 FflasLower = 122 63 }; 64 enum FFLAS_C_DIAG { 65 FflasNonUnit = 131, 66 FflasUnit = 132 67 }; 68 enum FFLAS_C_SIDE { 69 FflasLeft = 141, 70 FflasRight = 142 71 }; 72 73 #endif // __FFLASFFPACK_interfaces_libs_fflas_c_H 74 75 enum FFPACK_C_LU_TAG 76 { 77 FfpackSlabRecursive = 1, 78 FfpackTileRecursive = 2, 79 FfpackSingular = 3 80 }; 81 82 enum FFPACK_C_CHARPOLY_TAG 83 { 84 FfpackLUK=1, 85 FfpackKG=2, 86 FfpackHybrid=3, 87 FfpackKGFast=4, 88 FfpackDanilevski=5, 89 FfpackArithProg=6, 90 FfpackKGFastG=7 91 }; 92 93 enum FFPACK_C_MINPOLY_TAG 94 { 95 FfpackDense=1, 96 FfpackKGF=2 97 }; 98 99 100 101 /*****************/ 102 /* PERMUTATIONS */ 103 /*****************/ 104 105 106 void LAPACKPerm2MathPerm (size_t * MathP, const size_t * LapackP, 107 const size_t N); 108 109 void MathPerm2LAPACKPerm (size_t * LapackP, const size_t * MathP, 110 const size_t N); 111 112 void MatrixApplyS_modular_double (const double p, double * A, const size_t lda, const size_t width, 113 const size_t M2, 114 const size_t R1, const size_t R2, 115 const size_t R3, const size_t R4 116 , bool positive ); 117 118 void PermApplyS_double (double * A, const size_t lda, const size_t width, 119 const size_t M2, 120 const size_t R1, const size_t R2, 121 const size_t R3, const size_t R4); 122 123 124 void MatrixApplyT_modular_double (const double p, double * A, const size_t lda, const size_t width, 125 const size_t N2, 126 const size_t R1, const size_t R2, 127 const size_t R3, const size_t R4 128 , bool positive ); 129 130 void PermApplyT_double (double * A, const size_t lda, const size_t width, 131 const size_t N2, 132 const size_t R1, const size_t R2, 133 const size_t R3, const size_t R4); 134 135 136 void composePermutationsLLM (size_t * MathP, 137 const size_t * P1, 138 const size_t * P2, 139 const size_t R, const size_t N); 140 141 void composePermutationsLLL (size_t * P1, 142 const size_t * P2, 143 const size_t R, const size_t N); 144 145 void composePermutationsMLM (size_t * MathP1, 146 const size_t * P2, 147 const size_t R, const size_t N); 148 149 void cyclic_shift_mathPerm (size_t * P, const size_t s); 150 151 #if 0 152 template<typename Base_t> 153 void cyclic_shift_row_col(Base_t * A, size_t m, size_t n, size_t lda); 154 #endif 155 156 157 void cyclic_shift_row_modular_double(const double p, double * A, size_t m, size_t n, size_t lda 158 , bool positive ); 159 160 161 void cyclic_shift_col_modular_double(const double p, double * A, size_t m, size_t n, size_t lda 162 , bool positive ); 163 164 165 166 void 167 applyP_modular_double( const double p, 168 const enum FFLAS_C_SIDE Side, 169 const enum FFLAS_C_TRANSPOSE Trans, 170 const size_t M, const size_t ibeg, const size_t iend, 171 double * A, const size_t lda, const size_t * P 172 , bool positive ); 173 174 175 176 177 178 /* fgetrs, fgesv */ 179 180 void 181 fgetrsin_modular_double (const double p, 182 const enum FFLAS_C_SIDE Side, 183 const size_t M, const size_t N, const size_t R, 184 double * A, const size_t lda, 185 const size_t *P, const size_t *Q, 186 double * B, const size_t ldb, 187 int * info 188 , bool positive ); 189 190 191 double * 192 fgetrs_modular_double (const double p, 193 const enum FFLAS_C_SIDE Side, 194 const size_t M, const size_t N, const size_t NRHS, const size_t R, 195 double * A, const size_t lda, 196 const size_t *P, const size_t *Q, 197 double * X, const size_t ldx, 198 const double * B, const size_t ldb, 199 int * info 200 , bool positive ); 201 202 203 size_t 204 fgesvin_modular_double (const double p, 205 const enum FFLAS_C_SIDE Side, 206 const size_t M, const size_t N, 207 double * A, const size_t lda, 208 double * B, const size_t ldb, 209 int * info 210 , bool positive ); 211 212 213 size_t 214 fgesv_modular_double (const double p, 215 const enum FFLAS_C_SIDE Side, 216 const size_t M, const size_t N, const size_t NRHS, 217 double * A, const size_t lda, 218 double * X, const size_t ldx, 219 const double * B, const size_t ldb, 220 int * info); 221 222 /* ftrtr */ 223 224 225 void 226 ftrtri_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, const enum FFLAS_C_DIAG Diag, 227 const size_t N, double * A, const size_t lda 228 , bool positive ); 229 230 231 void trinv_left_modular_double( const double p, const size_t N, const double * L, const size_t ldl, 232 double * X, const size_t ldx 233 , bool positive ); 234 235 236 void 237 ftrtrm_modular_double (const double p, const enum FFLAS_C_DIAG diag, const size_t N, 238 double * A, const size_t lda 239 , bool positive ); 240 241 242 243 /* PLUQ */ 244 245 size_t 246 PLUQ_modular_double (const double p, const enum FFLAS_C_DIAG Diag, 247 const size_t M, const size_t N, 248 double * A, const size_t lda, 249 size_t*P, size_t *Q, bool positive); 250 251 252 size_t 253 LUdivine_modular_double (const double p, const enum FFLAS_C_DIAG Diag, const enum FFLAS_C_TRANSPOSE trans, 254 const size_t M, const size_t N, 255 double * A, const size_t lda, 256 size_t* P, size_t* Qt, 257 const enum FFPACK_C_LU_TAG LuTag, 258 const size_t cutoff 259 , bool positive ); 260 261 262 size_t 263 LUdivine_small_modular_double (const double p, const enum FFLAS_C_DIAG Diag, const enum FFLAS_C_TRANSPOSE trans, 264 const size_t M, const size_t N, 265 double * A, const size_t lda, 266 size_t* P, size_t* Q, 267 const enum FFPACK_C_LU_TAG LuTag 268 , bool positive ); 269 270 271 size_t 272 LUdivine_gauss_modular_double (const double p, const enum FFLAS_C_DIAG Diag, 273 const size_t M, const size_t N, 274 double * A, const size_t lda, 275 size_t* P, size_t* Q, 276 const enum FFPACK_C_LU_TAG LuTag 277 , bool positive ); 278 279 280 281 /*****************/ 282 /* ECHELON FORMS */ 283 /*****************/ 284 285 size_t 286 ColumnEchelonForm_modular_double (const double p, const size_t M, const size_t N, 287 double * A, const size_t lda, 288 size_t* P, size_t* Qt, bool transform, 289 const enum FFPACK_C_LU_TAG LuTag 290 , bool positive ); 291 292 293 size_t 294 RowEchelonForm_modular_double (const double p, const size_t M, const size_t N, 295 double * A, const size_t lda, 296 size_t* P, size_t* Qt, const bool transform, 297 const enum FFPACK_C_LU_TAG LuTag 298 , bool positive ); 299 300 size_t 301 ColumnEchelonForm_modular_float (const float p, const size_t M, const size_t N, 302 float * A, const size_t lda, 303 size_t* P, size_t* Qt, bool transform, 304 const enum FFPACK_C_LU_TAG LuTag 305 , bool positive ); 306 307 308 size_t 309 RowEchelonForm_modular_float (const float p, const size_t M, const size_t N, 310 float * A, const size_t lda, 311 size_t* P, size_t* Qt, const bool transform, 312 const enum FFPACK_C_LU_TAG LuTag 313 , bool positive ); 314 315 316 size_t 317 ColumnEchelonForm_modular_int32_t (const int32_t p, const size_t M, const size_t N, 318 int32_t * A, const size_t lda, 319 size_t* P, size_t* Qt, bool transform, 320 const enum FFPACK_C_LU_TAG LuTag 321 , bool positive ); 322 323 324 size_t 325 RowEchelonForm_modular_int32_t (const int32_t p, const size_t M, const size_t N, 326 int32_t * A, const size_t lda, 327 size_t* P, size_t* Qt, const bool transform, 328 const enum FFPACK_C_LU_TAG LuTag 329 , bool positive ); 330 331 332 size_t 333 ReducedColumnEchelonForm_modular_double (const double p, const size_t M, const size_t N, 334 double * A, const size_t lda, 335 size_t* P, size_t* Qt, const bool transform, 336 const enum FFPACK_C_LU_TAG LuTag 337 , bool positive ); 338 339 340 size_t 341 ReducedRowEchelonForm_modular_double (const double p, const size_t M, const size_t N, 342 double * A, const size_t lda, 343 size_t* P, size_t* Qt, const bool transform, 344 const enum FFPACK_C_LU_TAG LuTag 345 , bool positive ); 346 size_t 347 ReducedColumnEchelonForm_modular_float (const float p, const size_t M, const size_t N, 348 float * A, const size_t lda, 349 size_t* P, size_t* Qt, const bool transform, 350 const enum FFPACK_C_LU_TAG LuTag 351 , bool positive ); 352 353 354 size_t 355 ReducedRowEchelonForm_modular_float (const float p, const size_t M, const size_t N, 356 float * A, const size_t lda, 357 size_t* P, size_t* Qt, const bool transform, 358 const enum FFPACK_C_LU_TAG LuTag 359 , bool positive ); 360 361 size_t 362 ReducedColumnEchelonForm_modular_int32_t (const int32_t p, const size_t M, const size_t N, 363 int32_t * A, const size_t lda, 364 size_t* P, size_t* Qt, const bool transform, 365 const enum FFPACK_C_LU_TAG LuTag 366 , bool positive ); 367 368 369 size_t 370 ReducedRowEchelonForm_modular_int32_t (const int32_t p, const size_t M, const size_t N, 371 int32_t * A, const size_t lda, 372 size_t* P, size_t* Qt, const bool transform, 373 const enum FFPACK_C_LU_TAG LuTag 374 , bool positive ); 375 376 377 size_t 378 ReducedRowEchelonForm2_modular_double (const double p, const size_t M, const size_t N, 379 double * A, const size_t lda, 380 size_t* P, size_t* Qt, const bool transform 381 , bool positive ); 382 383 384 size_t 385 REF_modular_double (const double p, const size_t M, const size_t N, 386 double * A, const size_t lda, 387 const size_t colbeg, const size_t rowbeg, const size_t colsize, 388 size_t* Qt, size_t* P 389 , bool positive ); 390 391 392 /*****************/ 393 /* INVERSION */ 394 /*****************/ 395 396 397 double * 398 Invertin_modular_double (const double p, const size_t M, 399 double * A, const size_t lda, 400 int * nullity 401 , bool positive ); 402 403 404 double * 405 Invert_modular_double (const double p, const size_t M, 406 const double * A, const size_t lda, 407 double * X, const size_t ldx, 408 int* nullity 409 , bool positive ); 410 411 412 double * 413 Invert2_modular_double( const double p, const size_t M, 414 double * A, const size_t lda, 415 double * X, const size_t ldx, 416 int* nullity 417 , bool positive ); 418 419 /*****************************/ 420 /* CHARACTERISTIC POLYNOMIAL */ 421 /*****************************/ 422 423 424 #if 0 /* pas pour le moment */ 425 template <class Polynomial> 426 std::list<Polynomial>& 427 CharPoly( const double p, std::list<Polynomial>& charp, const size_t N, 428 double * A, const size_t lda, 429 const enum FFPACK_C_CHARPOLY_TAG CharpTag= FfpackArithProg); 430 431 template<class Polynomial> 432 Polynomial & mulpoly_modular_double(const double p, Polynomial &res, const Polynomial & P1, const Polynomial & P2); 433 434 template <class Polynomial> 435 Polynomial& 436 CharPoly_modular_double( const double p, Polynomial& charp, const size_t N, 437 double * A, const size_t lda, 438 const enum FFPACK_C_CHARPOLY_TAG CharpTag= FfpackArithProg); 439 440 441 442 template <class Polynomial> 443 std::list<Polynomial>& 444 CharpolyArithProg_modular_double (const double p, std::list<Polynomial>& frobeniusForm, 445 const size_t N, double * A, const size_t lda, const size_t c); 446 #endif 447 448 449 450 /**********************/ 451 /* MINIMAL POLYNOMIAL */ 452 /**********************/ 453 454 #if 0 /* pas pour le moment */ 455 template <class Polynomial> 456 Polynomial& 457 MinPoly_modular_double( const double p, Polynomial& minP, const size_t N, 458 const double * A, const size_t lda, 459 double * X, const size_t ldx, size_t* P, 460 const enum FFPACK_C_MINPOLY_TAG MinTag= FfpackDense, 461 const size_t kg_mc=0, const size_t kg_mb=0, const size_t kg_j=0 ); 462 #endif 463 464 465 /* Krylov Elim */ 466 467 468 size_t KrylovElim_modular_double( const double p, const size_t M, const size_t N, 469 double * A, const size_t lda, size_t*P, 470 size_t *Q, const size_t deg, size_t *iterates, size_t * inviterates, const size_t maxit,size_t virt 471 , bool positive ); 472 473 474 size_t SpecRankProfile_modular_double (const double p, const size_t M, const size_t N, 475 double * A, const size_t lda, const size_t deg, size_t *rankProfile 476 , bool positive ); 477 478 479 /********/ 480 /* RANK */ 481 /********/ 482 483 484 size_t 485 Rank_modular_double( const double p, const size_t M, const size_t N, 486 double * A, const size_t lda 487 , bool positive ) ; 488 489 /********/ 490 /* DET */ 491 /********/ 492 493 494 bool 495 IsSingular_modular_double( const double p, const size_t M, const size_t N, 496 double * A, const size_t lda 497 , bool positive ); 498 499 500 double 501 Det_modular_double( const double p, const size_t N, 502 double * A, const size_t lda 503 , bool positive ); 504 505 /*********/ 506 /* SOLVE */ 507 /*********/ 508 509 510 511 double * 512 Solve_modular_double( const double p, const size_t M, 513 double * A, const size_t lda, 514 double * x, const int incx, 515 const double * b, const int incb 516 , bool positive ); 517 518 519 520 void 521 solveLB_modular_double( const double p, const enum FFLAS_C_SIDE Side, 522 const size_t M, const size_t N, const size_t R, 523 double * L, const size_t ldl, 524 const size_t * Q, 525 double * B, const size_t ldb ); 526 527 528 void 529 solveLB2_modular_double( const double p, const enum FFLAS_C_SIDE Side, 530 const size_t M, const size_t N, const size_t R, 531 double * L, const size_t ldl, 532 const size_t * Q, 533 double * B, const size_t ldb 534 , bool positive ); 535 536 537 /*************/ 538 /* NULLSPACE */ 539 /*************/ 540 541 542 void RandomNullSpaceVector_modular_double (const double p, const enum FFLAS_C_SIDE Side, 543 const size_t M, const size_t N, 544 double * A, const size_t lda, 545 double * X, const size_t incX 546 , bool positive ); 547 548 549 size_t NullSpaceBasis_modular_double (const double p, const enum FFLAS_C_SIDE Side, 550 const size_t M, const size_t N, 551 double * A, const size_t lda, 552 double ** NS, size_t* ldn, 553 size_t * NSdim 554 , bool positive ); 555 556 /*****************/ 557 /* RANK PROFILES */ 558 /*****************/ 559 560 561 size_t RowRankProfile_modular_double (const double p, const size_t M, const size_t N, 562 double * A, const size_t lda, 563 size_t ** rkprofile, 564 const enum FFPACK_C_LU_TAG LuTag 565 , bool positive ); 566 567 568 569 size_t ColumnRankProfile_modular_double (const double p, const size_t M, const size_t N, 570 double * A, const size_t lda, 571 size_t ** rkprofile, 572 const enum FFPACK_C_LU_TAG LuTag 573 , bool positive ); 574 575 void RankProfileFromLU (const size_t* P, const size_t N, const size_t R, 576 size_t* rkprofile, const enum FFPACK_C_LU_TAG LuTag); 577 578 size_t LeadingSubmatrixRankProfiles (const size_t M, const size_t N, const size_t R, 579 const size_t LSm, const size_t LSn, 580 const size_t* P, const size_t* Q, 581 size_t* RRP, size_t* CRP); 582 583 584 size_t RowRankProfileSubmatrixIndices_modular_double (const double p, 585 const size_t M, const size_t N, 586 double * A, 587 const size_t lda, 588 size_t ** rowindices, 589 size_t ** colindices, 590 size_t * R 591 , bool positive ); 592 593 594 size_t ColRankProfileSubmatrixIndices_modular_double (const double p, 595 const size_t M, const size_t N, 596 double * A, 597 const size_t lda, 598 size_t** rowindices, 599 size_t** colindices, 600 size_t* R 601 , bool positive ); 602 603 604 size_t RowRankProfileSubmatrix_modular_double (const double p, 605 const size_t M, const size_t N, 606 double * A, 607 const size_t lda, 608 double ** X, size_t* R 609 , bool positive ); 610 611 612 size_t ColRankProfileSubmatrix_modular_double (const double p, const size_t M, const size_t N, 613 double * A, const size_t lda, 614 double ** X, size_t* R 615 , bool positive ); 616 617 /*********************************************/ 618 /* Accessors to Triangular and Echelon forms */ 619 /*********************************************/ 620 621 622 void 623 getTriangular_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 624 const enum FFLAS_C_DIAG diag, 625 const size_t M, const size_t N, const size_t R, 626 const double * A, const size_t lda, 627 double * T, const size_t ldt, 628 const bool OnlyNonZeroVectors 629 , bool positive ); 630 631 632 void 633 getTriangularin_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 634 const enum FFLAS_C_DIAG diag, 635 const size_t M, const size_t N, const size_t R, 636 double * A, const size_t lda 637 , bool positive ); 638 639 640 void 641 getEchelonForm_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 642 const enum FFLAS_C_DIAG diag, 643 const size_t M, const size_t N, const size_t R, const size_t* P, 644 const double * A, const size_t lda, 645 double * T, const size_t ldt, 646 const bool OnlyNonZeroVectors, 647 const enum FFPACK_C_LU_TAG LuTag 648 , bool positive ); 649 650 651 void 652 getEchelonFormin_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 653 const enum FFLAS_C_DIAG diag, 654 const size_t M, const size_t N, const size_t R, const size_t* P, 655 double * A, const size_t lda, 656 const enum FFPACK_C_LU_TAG LuTag 657 , bool positive ); 658 659 void 660 getEchelonTransform_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 661 const enum FFLAS_C_DIAG diag, 662 const size_t M, const size_t N, const size_t R, const size_t* P, const size_t* Q, 663 const double * A, const size_t lda, 664 double * T, const size_t ldt, 665 const enum FFPACK_C_LU_TAG LuTag 666 , bool positive ); 667 668 669 void 670 getReducedEchelonForm_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 671 const size_t M, const size_t N, const size_t R, const size_t* P, 672 const double * A, const size_t lda, 673 double * T, const size_t ldt, 674 const bool OnlyNonZeroVectors, 675 const enum FFPACK_C_LU_TAG LuTag 676 , bool positive ); 677 678 679 void 680 getReducedEchelonFormin_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 681 const size_t M, const size_t N, const size_t R, const size_t* P, 682 double * A, const size_t lda, 683 const enum FFPACK_C_LU_TAG LuTag 684 , bool positive ); 685 686 687 void 688 getReducedEchelonTransform_modular_double (const double p, const enum FFLAS_C_UPLO Uplo, 689 const size_t M, const size_t N, const size_t R, const size_t* P, const size_t* Q, 690 const double * A, const size_t lda, 691 double * T, const size_t ldt, 692 const enum FFPACK_C_LU_TAG LuTag 693 , bool positive ); 694 695 void 696 PLUQtoEchelonPermutation (const size_t N, const size_t R, const size_t * P, size_t * outPerm); 697 698 #ifdef __cplusplus 699 } 700 701 #endif 702 703 704 #endif // __FFLASFFPACK_interfaces_libs_ffpack_c_H 705 706 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 707 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 708