| | 686 | |
| | 687 | /** NullSpaceBasis |
| | 688 | * Computes a basis of the Left/Right nullspace of the matrix A |
| | 689 | * return the dimension of the nullspace. |
| | 690 | * |
| | 691 | * @param A: input matrix of dimension M x N, A is modified |
| | 692 | * @param NS: output matrix of dimension N x NSdim (allocated here) |
| | 693 | * @param NSdim: the dimension of the Nullspace (N-rank(A)) |
| | 694 | * |
| | 695 | */ |
| | 696 | template <class Field> |
| | 697 | static size_t NullSpaceBasis (const Field& F, const FFLAS_SIDE Side, |
| | 698 | const size_t M, const size_t N, |
| | 699 | typename Field::Element* A, const size_t lda, |
| | 700 | typename Field::Element*& NS, size_t& ldn, |
| | 701 | size_t& NSdim){ |
| | 702 | |
| | 703 | typename Field::Element one, zero,mone; |
| | 704 | F.init(one, 1.0); |
| | 705 | F.init(zero, 0.0); |
| | 706 | F.neg(mone, one); |
| | 707 | |
| | 708 | if (Side == FflasRight) { // Right NullSpace |
| | 709 | size_t* P = new size_t[N]; |
| | 710 | size_t* Qt = new size_t[M]; |
| | 711 | |
| | 712 | size_t R = LUdivine (F, FflasNonUnit, FflasNoTrans, M, N, A, lda, P, Qt); |
| | 713 | |
| | 714 | ldn = N-R; |
| | 715 | NSdim = ldn; |
| | 716 | NS = new typename Field::Element [N*ldn]; |
| | 717 | |
| | 718 | for (size_t i=0; i<R; ++i) |
| | 719 | fcopy (F, ldn, NS + i*ldn, 1, A + R + i*lda, 1); |
| | 720 | |
| | 721 | ftrsm (F, FflasLeft, FflasUpper, FflasNoTrans, FflasNonUnit, R, ldn, |
| | 722 | mone, A, lda, NS, ldn); |
| | 723 | |
| | 724 | for (size_t i=R; i<N; ++i){ |
| | 725 | for (size_t j=0; j < ldn; ++j) |
| | 726 | F.assign (*(NS+i*ldn+j), zero); |
| | 727 | F.assign (*(NS + i*ldn + i-R), one); |
| | 728 | } |
| | 729 | applyP (F, FflasLeft, FflasTrans, NSdim, 0, R, NS, ldn, P); |
| | 730 | delete [] P; |
| | 731 | delete [] Qt; |
| | 732 | return N-R; |
| | 733 | } else { // Left NullSpace |
| | 734 | size_t* P = new size_t[M]; |
| | 735 | size_t* Qt = new size_t[N]; |
| | 736 | |
| | 737 | size_t R = LUdivine (F, FflasNonUnit, FflasTrans, M, N, A, lda, P, Qt); |
| | 738 | |
| | 739 | write_field (F,cerr<<"LU="<<endl,A,M,N,lda); |
| | 740 | ldn = M; |
| | 741 | NSdim = M-R; |
| | 742 | NS = new typename Field::Element [NSdim*ldn]; |
| | 743 | for (size_t i=0; i<NSdim; ++i) |
| | 744 | fcopy (F, R, NS + i*ldn, 1, A + (R + i)*lda, 1); |
| | 745 | ftrsm (F, FflasRight, FflasLower, FflasNoTrans, FflasNonUnit, NSdim, R, |
| | 746 | mone, A, lda, NS, ldn); |
| | 747 | |
| | 748 | for (size_t i=0; i<NSdim; ++i){ |
| | 749 | for (size_t j=R; j < M; ++j) |
| | 750 | F.assign (*(NS+i*ldn+j), zero); |
| | 751 | F.assign (*(NS + i*ldn + i+R), one); |
| | 752 | } |
| | 753 | std::cerr<<"NSdim, M = "<<NSdim<<" "<<M<<std::endl; |
| | 754 | write_field (F,cerr<<"NS="<<endl,NS,NSdim,M,ldn); |
| | 755 | applyP (F, FflasRight, FflasNoTrans, NSdim, 0, R, NS, ldn, P); |
| | 756 | delete [] P; |
| | 757 | delete [] Qt; |
| | 758 | return N-R; |
| | 759 | } |
| | 760 | } |
| | 761 | |
