1/* ffpack_rankprofiles.inl 2 * Copyright (C) 2015 FFLAS-FFACK group 3 * 4 * Written by Clement Pernet <Clement.Pernet@imag.fr> 5 * 6 * ========LICENCE======== 7 * This file is part of the library FFLAS-FFPACK. 8 * 9 * FFLAS-FFPACK is free software: you can redistribute it and/or modify 10 * it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * This library is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with this library; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * ========LICENCE======== 23 *. 24 */ 25 26#ifndef __FFLASFFPACK_ffpack_rank_profiles_INL 27#define __FFLASFFPACK_ffpack_rank_profiles_INL 28 29namespace FFPACK{ 30 template <class Field> 31 inline size_t RowRankProfile (const Field& F, const size_t M, const size_t N, 32 typename Field::Element_ptr A, const size_t lda, 33 size_t* &rkprofile, const FFPACK_LU_TAG LuTag){ 34 FFLAS::ParSeqHelper::Sequential seqH; 35 return FFPACK::RowRankProfile (F, M, N, A, lda, rkprofile, LuTag, seqH); 36 } 37 38 template <class Field> 39 inline size_t pRowRankProfile (const Field& F, const size_t M, const size_t N, 40 typename Field::Element_ptr A, const size_t lda, 41 size_t* &rkprofile, size_t numthreads, const FFPACK_LU_TAG LuTag){ 42 size_t r; 43 PAR_BLOCK{ 44 size_t nt = numthreads ? numthreads : NUM_THREADS; 45 FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::Threads> parH(nt); 46 r = FFPACK::RowRankProfile (F, M, N, A, lda, rkprofile, LuTag, parH); 47 } 48 return r; 49 } 50 51 template <class Field, class PSHelper> 52 inline size_t RowRankProfile (const Field& F, const size_t M, const size_t N, 53 typename Field::Element_ptr A, const size_t lda, 54 size_t* &rkprofile, const FFPACK_LU_TAG LuTag, PSHelper& psH){ 55 56 57 size_t *P = FFLAS::fflas_new<size_t>((LuTag==FfpackSlabRecursive)?N:M); 58 size_t *Q = FFLAS::fflas_new<size_t>((LuTag==FfpackSlabRecursive)?M:N); 59 size_t R; 60 61 if (LuTag == FfpackSlabRecursive){ 62 R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasNoTrans, M, N, A, lda, P, Q); 63 std::swap(P,Q); 64 } else 65 R = PLUQ (F, FFLAS::FflasNonUnit, M, N, A, lda, P, Q, psH); 66 67 rkprofile = FFLAS::fflas_new<size_t> (R); 68 69 RankProfileFromLU (P, M, R, rkprofile, LuTag); 70 71 FFLAS::fflas_delete (Q); 72 FFLAS::fflas_delete (P); 73 return R; 74 } 75 76 template <class Field> 77 inline size_t ColumnRankProfile (const Field& F, const size_t M, const size_t N, 78 typename Field::Element_ptr A, const size_t lda, 79 size_t* &rkprofile, const FFPACK_LU_TAG LuTag){ 80 FFLAS::ParSeqHelper::Sequential seqH; 81 return FFPACK::ColumnRankProfile (F, M, N, A, lda, rkprofile, LuTag, seqH); 82 } 83 84 template <class Field> 85 inline size_t pColumnRankProfile (const Field& F, const size_t M, const size_t N, 86 typename Field::Element_ptr A, const size_t lda, 87 size_t* &rkprofile, size_t numthreads, const FFPACK_LU_TAG LuTag){ 88 size_t r; 89 PAR_BLOCK{ 90 size_t nt = numthreads ? numthreads : NUM_THREADS; 91 FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Recursive,FFLAS::StrategyParameter::Threads> parH(nt); 92 r = FFPACK::ColumnRankProfile (F, M, N, A, lda, rkprofile, LuTag, parH); 93 } 94 return r; 95 } 96 97 template <class Field, class PSHelper> 98 inline size_t ColumnRankProfile (const Field& F, const size_t M, const size_t N, 99 typename Field::Element_ptr A, const size_t lda, 100 size_t* &rkprofile, const FFPACK_LU_TAG LuTag, PSHelper& psH){ 101 102 size_t *P = FFLAS::fflas_new<size_t>(M); 103 size_t *Q = FFLAS::fflas_new<size_t>(N); 104 size_t R; 105 106 if (LuTag == FfpackSlabRecursive){ 107 R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasTrans, M, N, A, lda, P, Q); 108 } else 109 R = PLUQ (F, FFLAS::FflasNonUnit, M, N, A, lda, P, Q, psH); 110 111 rkprofile = FFLAS::fflas_new<size_t> (R); 112 113 RankProfileFromLU (Q, N, R, rkprofile, LuTag); 114 115 FFLAS::fflas_delete (P); 116 FFLAS::fflas_delete (Q); 117 return R; 118 } 119 120 121 inline void RankProfileFromLU (const size_t* Q, const size_t N, const size_t R, 122 size_t* rkprofile, const FFPACK_LU_TAG LuTag){ 123 124 if (LuTag == FfpackSlabRecursive) 125 std::copy(Q, Q+R, rkprofile); 126 else { 127 size_t * RP = FFLAS::fflas_new<size_t>(N); 128 for (size_t i=0;i < N; ++i) 129 RP [i] = i; 130 for (size_t i=0; i<N; ++i) 131 if (Q[i] != i) 132 std::swap (RP [i], RP [Q [i]]); 133 134 std::copy(RP, RP+R, rkprofile); 135 std::sort (rkprofile, rkprofile + R); 136 FFLAS::fflas_delete(RP); 137 } 138 } 139 140 inline size_t LeadingSubmatrixRankProfiles (const size_t M, const size_t N, const size_t R, 141 const size_t LSm, const size_t LSn, 142 const size_t* P, const size_t* Q, 143 size_t* RRP, size_t* CRP){ 144 size_t LSr=0; // rank of the LSm x LSn leading submatrix 145 146 size_t* MathP = FFLAS::fflas_new<size_t>(M); 147 size_t* MathQ = FFLAS::fflas_new<size_t>(N); 148 149 LAPACKPerm2MathPerm (MathP, P, M); 150 LAPACKPerm2MathPerm (MathQ, Q, N); 151 for (size_t i = 0; i < R; i++) 152 if (MathP[i] < LSm && MathQ[i] < LSn){ 153 RRP [LSr] = MathP[i]; 154 CRP [LSr] = MathQ[i]; 155 LSr++; 156 } 157 std::sort (RRP, RRP+LSr); 158 std::sort (CRP, CRP+LSr); 159 FFLAS::fflas_delete(MathP); 160 FFLAS::fflas_delete(MathQ); 161 return LSr; 162 163 } 164 165 166 template <class Field> 167 size_t RowRankProfileSubmatrixIndices (const Field& F, 168 const size_t M, const size_t N, 169 typename Field::Element_ptr A, 170 const size_t lda, 171 size_t*& rowindices, 172 size_t*& colindices, 173 size_t& R) 174 { 175 size_t *P = FFLAS::fflas_new<size_t>(N); 176 size_t *Q = FFLAS::fflas_new<size_t>(M); 177 178 R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasNoTrans, M, N, A, lda, P, Q); 179 rowindices = FFLAS::fflas_new<size_t>(M); 180 colindices = FFLAS::fflas_new<size_t>(N); 181 for (size_t i=0; i<R; ++i){ 182 rowindices [i] = Q [i]; 183 } 184 for (size_t i=0; i<N; ++i) 185 colindices [i] = i; 186 size_t tmp; 187 for (size_t i=0; i<R; ++i){ 188 if (i != P[i]){ 189 tmp = colindices[i]; 190 colindices[i] = colindices[P[i]]; 191 colindices[P[i]] = tmp; 192 } 193 } 194 195 FFLAS::fflas_delete( P); 196 FFLAS::fflas_delete( Q); 197 198 return R; 199 } 200 201 template <class Field> 202 size_t ColRankProfileSubmatrixIndices (const Field& F, 203 const size_t M, const size_t N, 204 typename Field::Element_ptr A, 205 const size_t lda, 206 size_t*& rowindices, 207 size_t*& colindices, 208 size_t& R) 209 { 210 size_t *P = FFLAS::fflas_new<size_t>(M); 211 size_t *Q = FFLAS::fflas_new<size_t>(N); 212 213 R = LUdivine (F, FFLAS::FflasNonUnit, FFLAS::FflasTrans, M, N, A, lda, P, Q); 214 rowindices = FFLAS::fflas_new<size_t>(M); 215 colindices = FFLAS::fflas_new<size_t>(N); 216 for (size_t i=0; i<R; ++i) 217 colindices [i] = Q [i]; 218 219 for (size_t i=0; i<N; ++i) 220 rowindices [i] = i; 221 222 size_t tmp; 223 for (size_t i=0; i<R; ++i){ 224 if (i != P[i]){ 225 tmp = rowindices[i]; 226 rowindices[i] = rowindices[P[i]]; 227 rowindices[P[i]] = tmp; 228 } 229 } 230 FFLAS::fflas_delete( P); 231 FFLAS::fflas_delete( Q); 232 233 return R; 234 } 235 236 template <class Field> 237 size_t RowRankProfileSubmatrix (const Field& F, 238 const size_t M, const size_t N, 239 typename Field::Element_ptr A, 240 const size_t lda, 241 typename Field::Element_ptr& X, size_t& R) 242 { 243 244 size_t * rowindices, * colindices; 245 246 typename Field::Element_ptr A2 = FFLAS::fflas_new (F, M, N) ; 247 FFLAS::fassign(F,M,N,A,lda,A2,N); 248 249 RowRankProfileSubmatrixIndices (F, M, N, A2, N, rowindices, colindices, R); 250 251 X = FFLAS::fflas_new (F, R, R); 252 for (size_t i=0; i<R; ++i) 253 for (size_t j=0; j<R; ++j) 254 F.assign (*(X + i*R + j), *(A + rowindices[i]*lda + colindices[j])); 255 FFLAS::fflas_delete (A2); 256 FFLAS::fflas_delete( rowindices); 257 FFLAS::fflas_delete( colindices); 258 return R; 259 } 260 261 template <class Field> 262 size_t ColRankProfileSubmatrix (const Field& F, const size_t M, const size_t N, 263 typename Field::Element_ptr A, const size_t lda, 264 typename Field::Element_ptr& X, size_t& R) 265 { 266 267 size_t * rowindices, * colindices; 268 269 typename Field::Element_ptr A2 = FFLAS::fflas_new (F, M, N); 270 FFLAS::fassign(F,M,N,A,lda,A2,N); 271 272 ColRankProfileSubmatrixIndices (F, M, N, A2, N, rowindices, colindices, R); 273 274 X = FFLAS::fflas_new (F, R, R); 275 for (size_t i=0; i<R; ++i) 276 for (size_t j=0; j<R; ++j) 277 F.assign (*(X + i*R + j), *(A + rowindices[i]*lda + colindices[j])); 278 FFLAS::fflas_delete (A2); 279 FFLAS::fflas_delete( colindices); 280 FFLAS::fflas_delete( rowindices); 281 return R; 282 } 283 284 template <class Field> 285 typename Field::Element_ptr 286 LQUPtoInverseOfFullRankMinor( const Field& F, const size_t rank, 287 typename Field::Element_ptr A_factors, const size_t lda, 288 const size_t* QtPointer, 289 typename Field::Element_ptr X, const size_t ldx) 290 { 291 292 // upper entries are okay, just need to move up bottom ones 293 const size_t* srcRow = QtPointer; 294 for (size_t row=0; row<rank; row++, srcRow++) 295 if (*srcRow != row) { 296 typename Field::Element_ptr oldRow = A_factors + (*srcRow) * lda; 297 typename Field::Element_ptr newRow = A_factors + row * lda; 298 for (size_t col=0; col<row; col++, oldRow++, newRow++) 299 F.assign(*newRow, *oldRow); 300 } 301 302 // X <- (Qt.L.Q)^(-1) 303 //invL( F, rank, A_factors, lda, X, ldx); 304 ftrtri (F, FFLAS::FflasLower, FFLAS::FflasUnit, rank, A_factors, lda); 305 FFLAS::fassign(F,rank,rank,A_factors,lda,X,ldx); 306 307 // X = U^-1.X 308 ftrsm( F, FFLAS::FflasLeft, FFLAS::FflasUpper, FFLAS::FflasNoTrans, 309 FFLAS::FflasNonUnit, rank, rank, F.one, A_factors, lda, X, ldx); 310 311 return X; 312 313 } 314 315} // namespace FFPACK 316 317#endif // __FFLASFFPACK_ffpack_rank_profiles_INL 318/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 319// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 320