1 /* 2 * Copyright (C) FFLAS-FFPACK 3 * Written by Brice Boyer (briceboyer) <boyer.brice@gmail.com> 4 * This file is Free Software and part of FFLAS-FFPACK. 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 /* @file utils/fflas_randommatrix.h 27 * @ingroup tests 28 * @brief Utilities to create matrices with prescribed shapes, properties,... 29 * To be used in benchmarks/tests 30 */ 31 32 #ifndef __FFLASFFPACK_randommatrix_H 33 #define __FFLASFFPACK_randommatrix_H 34 35 #include "fflas-ffpack/fflas-ffpack-config.h" 36 #include "fflas-ffpack/utils/debug.h" 37 #include "fflas-ffpack/fflas/fflas.h" 38 #include <givaro/givinteger.h> 39 #include <givaro/givintprime.h> 40 #include <givaro/givranditer.h> 41 42 namespace FFPACK { 43 44 /** @brief Random non-zero Matrix. 45 * Creates a \c m x \c n matrix with random entries, and at least one of them is non zero. 46 * @param F field 47 * @param m number of rows in \p A 48 * @param n number of cols in \p A 49 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 50 * @param lda leading dimension of \p A 51 * @param G a random iterator 52 * @return \c A. 53 */ 54 template<class Field, class RandIter> 55 inline typename Field::Element_ptr NonZeroRandomMatrix(const Field & F,size_t m,size_t n,typename Field::Element_ptr A,size_t lda,RandIter & G)56 NonZeroRandomMatrix(const Field & F, size_t m, size_t n, typename Field::Element_ptr A, size_t lda, RandIter& G) { 57 bool ok=false; 58 while (!ok) 59 for (size_t i=0 ; i<m ; ++i) 60 for (size_t j= 0; j<n ;++j) 61 if (!F.isZero(G.random (A[i*lda+j]))) 62 ok = true; 63 return A; 64 } 65 66 /** @brief Random non-zero Matrix. 67 * Creates a \c m x \c n matrix with random entries, and at least one of them is non zero. 68 * @param F field 69 * @param m number of rows in \p A 70 * @param n number of cols in \p A 71 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 72 * @param lda leading dimension of \p A 73 * @return \c A. 74 */ 75 template<class Field, class RandIter> 76 inline typename Field::Element_ptr NonZeroRandomMatrix(const Field & F,size_t m,size_t n,typename Field::Element_ptr A,size_t lda)77 NonZeroRandomMatrix(const Field & F, size_t m, size_t n, 78 typename Field::Element_ptr A, size_t lda) { 79 typename Field::RandIter G(F); 80 return NonZeroRandomMatrix(F, m, n, A, lda, G); 81 } 82 83 /** @brief Random Matrix. 84 * Creates a \c m x \c n matrix with random entries. 85 * @param F field 86 * @param m number of rows in \p A 87 * @param n number of cols in \p A 88 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 89 * @param lda leading dimension of \p A 90 * @param G a random iterator 91 * @return \c A. 92 */ 93 template<class Field, class RandIter> 94 inline typename Field::Element_ptr RandomMatrix(const Field & F,size_t m,size_t n,typename Field::Element_ptr A,size_t lda,RandIter & G)95 RandomMatrix(const Field & F, size_t m, size_t n, typename Field::Element_ptr A, size_t lda, RandIter& G) { 96 for (size_t i=0 ; i<m ; ++i) 97 for (size_t j= 0; j<n ;++j) 98 G.random (A[i*lda+j]); 99 return A; 100 } 101 102 /** @brief Random Matrix. 103 * Creates a \c m x \c n matrix with random entries. 104 * @param F field 105 * @param m number of rows in \p A 106 * @param n number of cols in \p A 107 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 108 * @param lda leading dimension of \p A 109 * @return \c A. 110 */ 111 template<class Field> 112 inline typename Field::Element_ptr RandomMatrix(const Field & F,size_t m,size_t n,typename Field::Element_ptr A,size_t lda)113 RandomMatrix(const Field & F, size_t m, size_t n, typename Field::Element_ptr A, size_t lda) { 114 typename Field::RandIter G(F); 115 return RandomMatrix (F, m, n, A, lda, G); 116 } 117 118 /** @brief Random Triangular Matrix. 119 * Creates a \c m x \c n triangular matrix with random entries. The \c UpLo parameter defines wether it is upper or lower triangular. 120 * @param F field 121 * @param m number of rows in \p A 122 * @param n number of cols in \p A 123 * @param UpLo whether \c A is upper or lower triangular 124 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 125 * @param lda leading dimension of \p A 126 * @param G a random iterator 127 * @return \c A. 128 */ 129 template<class Field, class RandIter> 130 inline typename Field::Element_ptr RandomTriangularMatrix(const Field & F,size_t m,size_t n,const FFLAS::FFLAS_UPLO UpLo,const FFLAS::FFLAS_DIAG Diag,bool nonsingular,typename Field::Element_ptr A,size_t lda,RandIter & G)131 RandomTriangularMatrix (const Field & F, size_t m, size_t n, 132 const FFLAS::FFLAS_UPLO UpLo, const FFLAS::FFLAS_DIAG Diag, bool nonsingular, 133 typename Field::Element_ptr A, size_t lda, RandIter& G) { 134 135 if (UpLo == FFLAS::FflasUpper){ 136 for (size_t i=0 ; i<m ; ++i){ 137 FFLAS::fzero(F, std::min(i,n), A + i*lda, 1); 138 for (size_t j= i; j<n ;++j) 139 G.random (A[i*lda+j]); 140 } 141 } else { // FflasLower 142 for (size_t i=0 ; i<m ; ++i){ 143 for (size_t j=0; j<=i ;++j) 144 G.random (A[i*lda+j]); 145 FFLAS::fzero (F, n-1-std::min(i,n-1), A + i*lda + i+1, 1); 146 } 147 } 148 if (Diag == FFLAS::FflasUnit){ 149 for (size_t i=0; i< std::min(m,n); ++i) 150 F.assign(A[i*(lda+1)], F.one); 151 } else { // NonUnit 152 if (nonsingular){ 153 Givaro::GeneralRingNonZeroRandIter<Field,RandIter> nzG (G); 154 for (size_t i=0; i< std::min(m,n); ++i) 155 nzG.random (A[i*(lda+1)]); 156 } 157 } 158 return A; 159 } 160 /** @brief Random Triangular Matrix. 161 * Creates a \c m x \c n triangular matrix with random entries. The \c UpLo parameter defines wether it is upper or lower triangular. 162 * @param F field 163 * @param m number of rows in \p A 164 * @param n number of cols in \p A 165 * @param UpLo whether \c A is upper or lower triangular 166 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 167 * @param lda leading dimension of \p A 168 * @return \c A. 169 */ 170 template<class Field> 171 inline typename Field::Element_ptr RandomTriangularMatrix(const Field & F,size_t m,size_t n,const FFLAS::FFLAS_UPLO UpLo,const FFLAS::FFLAS_DIAG Diag,bool nonsingular,typename Field::Element_ptr A,size_t lda)172 RandomTriangularMatrix (const Field & F, size_t m, size_t n, 173 const FFLAS::FFLAS_UPLO UpLo, const FFLAS::FFLAS_DIAG Diag, bool nonsingular, 174 typename Field::Element_ptr A, size_t lda) { 175 typename Field::RandIter G(F); 176 return RandomTriangularMatrix (F, m, n, UpLo, Diag, nonsingular, A, lda, G); 177 } 178 179 /* Random integer in range. 180 * @param a min bound 181 * @param b max bound 182 * @return a random integer in [a,b[ */ RandInt(size_t a,size_t b)183 inline size_t RandInt(size_t a, size_t b) 184 { 185 size_t x = a ; 186 x += (size_t)rand()%(b-a); 187 FFLASFFPACK_check(x<b && x>=a); 188 return x ; 189 } 190 /** @brief Random Symmetric Matrix. 191 * Creates a \c m x \c n triangular matrix with random entries. The \c UpLo parameter defines wether it is upper or lower triangular. 192 * @param F field 193 * @param n order of \p A 194 * @param [out] A the matrix (preallocated to at least \c n x \c lda field elements) 195 * @param lda leading dimension of \p A 196 * @param G a random iterator 197 * @return \c A. 198 */ 199 template<class Field, class RandIter> 200 inline typename Field::Element_ptr RandomSymmetricMatrix(const Field & F,size_t n,bool nonsingular,typename Field::Element_ptr A,size_t lda,RandIter & G)201 RandomSymmetricMatrix (const Field & F,size_t n, bool nonsingular, 202 typename Field::Element_ptr A, size_t lda, RandIter& G) { 203 RandomTriangularMatrix (F, n, n, FFLAS::FflasUpper, FFLAS::FflasNonUnit, nonsingular, A, lda, G); 204 for (size_t i=0; i<n; i++){ 205 typename Field::Element piv = A[i*(lda+1)]; 206 if (!F.isZero(piv)){ 207 typename Field::Element inv; 208 F.init(inv); 209 F.inv(inv, A[i*(lda+1)]); 210 FFLAS::fscal(F, n-i-1, inv, A+i*(lda+1)+1, 1, A+i*(lda+1)+lda, lda); 211 } 212 } 213 ftrtrm (F, FFLAS::FflasRight, FFLAS::FflasNonUnit, n, A, lda); 214 return A; 215 } 216 } // FFPACK 217 218 #include "fflas-ffpack/ffpack/ffpack.h" 219 220 namespace FFPACK{ 221 /** @brief Random Matrix with prescribed rank. 222 * Creates an \c m x \c n matrix with random entries and rank \c r. 223 * @param F field 224 * @param m number of rows in \p A 225 * @param n number of cols in \p A 226 * @param r rank of the matrix to build 227 * @param A the matrix (preallocated to at least \c m x \c lda field elements) 228 * @param lda leading dimension of \p A 229 * @param G a random iterator 230 * @return \c A. 231 */ 232 template<class Field, class RandIter> 233 inline typename Field::Element_ptr RandomMatrixWithRank(const Field & F,size_t m,size_t n,size_t r,typename Field::Element_ptr A,size_t lda,RandIter & G)234 RandomMatrixWithRank (const Field & F, size_t m, size_t n, size_t r, 235 typename Field::Element_ptr A, size_t lda, RandIter& G){ 236 FFLASFFPACK_check(r <= std::min(m,n)); 237 FFLASFFPACK_check(n <= lda); 238 typedef typename Field::Element_ptr Element_ptr; 239 240 Givaro::GeneralRingNonZeroRandIter<Field,RandIter> nzG (G); 241 242 size_t * P = FFLAS::fflas_new<size_t>(n); 243 size_t * Q = FFLAS::fflas_new<size_t>(m); 244 for (size_t i = 0 ; i < m ; ++i ) Q[i] = 0; 245 for (size_t i = 0 ; i < n ; ++i ) P[i] = 0; 246 247 Element_ptr U = FFLAS::fflas_new(F,m,n); 248 Element_ptr L = FFLAS::fflas_new(F,m,m); 249 /* Create L, lower invertible */ 250 RandomTriangularMatrix (F, m, m, FFLAS::FflasLower, FFLAS::FflasNonUnit, true, L, m, G); 251 /* Create U, upper or rank r */ 252 RandomTriangularMatrix (F, m, n, FFLAS::FflasUpper, FFLAS::FflasNonUnit, true, U, n, G); 253 254 /* Create a random P,Q */ 255 for (size_t i = 0 ; i < n ; ++i) 256 P[i] = i + RandInt(0U,n-i); 257 for (size_t i = 0 ; i < m ; ++i) 258 Q[i] = i + RandInt(0U,m-i); 259 260 /* compute product */ 261 262 FFPACK::applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, 263 m,0,(int)n, U, n, P); 264 FFPACK::applyP (F, FFLAS::FflasLeft, FFLAS::FflasNoTrans, 265 m,0,(int)m, L, m, Q); 266 FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, 267 m, n, m, F.one, L, m, U, n, F.zero, A, lda); 268 // @todo compute LU with ftrtr 269 270 FFLAS::fflas_delete(P); 271 FFLAS::fflas_delete(L); 272 FFLAS::fflas_delete(U); 273 FFLAS::fflas_delete(Q); 274 275 return A; 276 } 277 278 /** @brief Random Matrix with prescribed rank. 279 * Creates an \c m x \c n matrix with random entries and rank \c r. 280 * @param F field 281 * @param m number of rows in \p A 282 * @param n number of cols in \p A 283 * @param r rank of the matrix to build 284 * @param [out] A the matrix (preallocated to at least \c m x \c lda field elements) 285 * @param lda leading dimension of \p A 286 * @return \c A. 287 */ 288 template<class Field> 289 inline typename Field::Element_ptr RandomMatrixWithRank(const Field & F,size_t m,size_t n,size_t r,typename Field::Element_ptr A,size_t lda)290 RandomMatrixWithRank (const Field & F, size_t m, size_t n, size_t r, 291 typename Field::Element_ptr A, size_t lda){ 292 typename Field::RandIter G(F); 293 return RandomMatrixWithRank(F, m, n, r, A, lda, G); 294 } 295 296 /** @brief Pick uniformly at random a sequence of \c R distinct elements from the set \f$ \{0,\dots, N-1\}\f$ using Knuth's shuffle. 297 * @param N the cardinality of the sampling set 298 * @param R the number of elements to sample 299 * @param [out] P the output sequence (pre-allocated to at least R indices) 300 */ RandomIndexSubset(size_t N,size_t R,size_t * P)301 inline size_t * RandomIndexSubset (size_t N, size_t R, size_t* P){ 302 size_t * Q = FFLAS::fflas_new<size_t>(N); 303 for (size_t i=0; i<N; ++i) 304 Q[i] = i; 305 for (size_t i=0; i<R; ++i){ 306 size_t j = RandInt(i,N); 307 P[i] = Q[j]; 308 Q[j] = Q[i]; 309 } 310 FFLAS::fflas_delete(Q); 311 return P; 312 } 313 314 /** @brief Pick uniformly at random a permutation of size \c N stored in LAPACK format using Knuth's shuffle. 315 * @param N the length of the permutation 316 * @param [out] P the output permutation (pre-allocated to at least N indices) 317 */ RandomPermutation(size_t N,size_t * P)318 inline size_t * RandomPermutation (size_t N, size_t* P){ 319 for (size_t i = 0 ; i < N ; ++i) 320 P[i] = i + RandInt(0U,N-i); 321 return P; 322 } 323 324 /** @brief Pick uniformly at random an R-subpermutation of dimension \c M x \c N : a matrix with only R non-zeros equal to one, in a random rook placement. 325 * @param M row dimension 326 * @param N column dimension 327 * @param [out] rows the row position of each non zero element (pre-allocated) 328 * @param [out] cols the column position of each non zero element (pre-allocated) 329 */ RandomRankProfileMatrix(size_t M,size_t N,size_t R,size_t * rows,size_t * cols)330 inline void RandomRankProfileMatrix (size_t M, size_t N, size_t R, size_t* rows, size_t* cols){ 331 RandomIndexSubset (M, R, rows); 332 RandomIndexSubset (N, R, cols); 333 } 334 swapval(size_t k,size_t N,size_t * P,size_t val)335 inline void swapval(size_t k, size_t N, size_t * P, size_t val){ 336 size_t i = k; 337 int found =-1; 338 do { 339 if (P[i] == val) 340 found = i; 341 i++; 342 } while(found<0); 343 P[found] = P[k]; 344 } 345 /** @brief Pick uniformly at random a symmetric R-subpermutation of dimension \c N x \c N : a symmetric matrix with only R non-zeros, all equal to one, in a random rook placement. 346 * @param N matrix order 347 * @param [out] rows the row position of each non zero element (pre-allocated) 348 * @param [out] cols the column position of each non zero element (pre-allocated) 349 */ RandomSymmetricRankProfileMatrix(size_t N,size_t R,size_t * rows,size_t * cols)350 inline void RandomSymmetricRankProfileMatrix (size_t N, size_t R, size_t* rows, size_t* cols){ 351 352 size_t * rr = FFLAS::fflas_new<size_t>(N); 353 size_t * cc = FFLAS::fflas_new<size_t>(N); 354 for (size_t i=0; i<N; ++i) 355 rr[i] = cc[i] = i; 356 for (size_t k=0; k<R; k+=2){ 357 size_t i = RandInt(k,N); 358 size_t j = RandInt(k,N); 359 cols[k] = cc[j]; 360 cc[j] = cc[k]; 361 rows[k] = rr[i]; 362 rr[i] = rr[k]; 363 if (rows[k] != cols[k] && k < R-1){ 364 // adding the symmetric element 365 rows[k+1] = cols[k]; 366 cols[k+1] = rows[k]; 367 swapval(k+1,N,rr,cols[k]); 368 swapval(k+1,N,cc,rows[k]); 369 } else { 370 // we need to add a diagonal pivot since 371 // - either k==R-1 and there is only one pivot left to be added 372 // - or we just added a diagonal pivot. We need to pick another one so 373 // that they appear with the same probability 2/N^2 as off-diagonal pivots 374 if (k<R-1) k++; // 375 size_t l, co; 376 int found =-1; 377 do{ 378 l = RandInt(k,N); 379 co = cc[l]; 380 for (size_t m=k; m<N; m++) 381 if (rr[m] == co) // l is valid as row co still available 382 found = m; 383 // TODO: Write a variant for when k < N/2 384 } while(found<0); 385 cols[k] = co; 386 cc[l] = cc[k]; 387 rows[k] = co; 388 rr[found] = rr[k]; 389 if (k<R) k--; // 390 } 391 } 392 FFLAS::fflas_delete(rr,cc); 393 } 394 395 /** @brief Random Matrix with prescribed rank and rank profile matrix 396 * Creates an \c m x \c n matrix with random entries and rank \c r. 397 * @param F field 398 * @param m number of rows in \p A 399 * @param n number of cols in \p A 400 * @param r rank of the matrix to build 401 * @param A the matrix (preallocated to at least \c m x \c lda field elements) 402 * @param lda leading dimension of \p A 403 * @param RRP the R dimensional array with row positions of the rank profile matrix' pivots 404 * @param CRP the R dimensional array with column positions of the rank profile matrix' pivots 405 * @param G a random iterator 406 * @return \c A. 407 */ 408 template<class Field,class RandIter> 409 inline typename Field::Element_ptr RandomMatrixWithRankandRPM(const Field & F,size_t M,size_t N,size_t R,typename Field::Element_ptr A,size_t lda,const size_t * RRP,const size_t * CRP,RandIter & G)410 RandomMatrixWithRankandRPM (const Field& F, size_t M, size_t N, size_t R, 411 typename Field::Element_ptr A, size_t lda, 412 const size_t * RRP, const size_t * CRP, RandIter& G){ 413 414 Givaro::GeneralRingNonZeroRandIter<Field,RandIter> nzG(G); 415 416 typename Field::Element_ptr L= FFLAS::fflas_new(F,M,N); 417 418 FFLAS::fzero(F, M, N, L, N); 419 // Disabling the parallel loop, as there is no way to declare G as SHARED in paladin 420 //FFLAS::ParSeqHelper::Parallel<FFLAS::CuttingStrategy::Block,FFLAS::StrategyParameter::Threads> H; 421 //SYNCH_GROUP (FOR1D(k, R, H, 422 for (size_t k=0; k<R; ++k){ 423 size_t i = RRP[k]; 424 size_t j = CRP[k]; 425 nzG.random (L [i*N+j]); 426 for (size_t l=i+1; l < M; ++l) 427 G.random (L [l*N+j]); 428 } 429 //)); 430 431 typename Field::Element_ptr U= FFLAS::fflas_new(F,N,N); 432 RandomTriangularMatrix (F, N, N, FFLAS::FflasUpper, FFLAS::FflasNonUnit, true, U, N, G); 433 434 // auto sp=SPLITTER(); //CP: broken with Modular<Integer>. Need to reorganize the helper behaviour with ParSeq and ModeTraits 435 auto sp=NOSPLIT(); 436 FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, M, N, N, 437 F.one, L, N, U, N, F.zero, A, lda, sp); 438 439 FFLAS::fflas_delete(L); 440 FFLAS::fflas_delete(U); 441 return A; 442 } 443 444 /** @brief Random Matrix with prescribed rank and rank profile matrix 445 * Creates an \c m x \c n matrix with random entries and rank \c r. 446 * @param F field 447 * @param m number of rows in \p A 448 * @param n number of cols in \p A 449 * @param r rank of the matrix to build 450 * @param A the matrix (preallocated to at least \c m x \c lda field elements) 451 * @param lda leading dimension of \p A 452 * @param RRP the R dimensional array with row positions of the rank profile matrix' pivots 453 * @param CRP the R dimensional array with column positions of the rank profile matrix' pivots 454 * @return \c A. 455 */ 456 template<class Field> 457 inline typename Field::Element_ptr RandomMatrixWithRankandRPM(const Field & F,size_t M,size_t N,size_t R,typename Field::Element_ptr A,size_t lda,const size_t * RRP,const size_t * CRP)458 RandomMatrixWithRankandRPM (const Field& F, size_t M, size_t N, size_t R, 459 typename Field::Element_ptr A, size_t lda, 460 const size_t * RRP, const size_t * CRP){ 461 typename Field::RandIter G(F); 462 return RandomMatrixWithRankandRPM (F, M, N, R, A, lda, RRP, CRP, G); 463 } 464 465 /** @brief Random Symmetric Matrix with prescribed rank and rank profile matrix 466 * Creates an \c n x \c n symmetric matrix with random entries and rank \c r. 467 * @param F field 468 * @param n order of \p A 469 * @param r rank of \p A 470 * @param A the matrix (preallocated to at least \c n x \c lda field elements) 471 * @param lda leading dimension of \p A 472 * @param RRP the R dimensional array with row positions of the rank profile matrix' pivots 473 * @param CRP the R dimensional array with column positions of the rank profile matrix' pivots 474 * @param G a random iterator 475 * @return \c A. 476 */ 477 template<class Field,class RandIter> 478 inline typename Field::Element_ptr RandomSymmetricMatrixWithRankandRPM(const Field & F,size_t N,size_t R,typename Field::Element_ptr A,size_t lda,const size_t * RRP,const size_t * CRP,RandIter & G)479 RandomSymmetricMatrixWithRankandRPM (const Field& F, size_t N, size_t R, 480 typename Field::Element_ptr A, size_t lda, 481 const size_t * RRP, const size_t * CRP, RandIter& G){ 482 483 typename Field::Element_ptr U= FFLAS::fflas_new(F,N,N); 484 typename Field::Element_ptr L= FFLAS::fflas_new(F,N,N); 485 // U <- $ 486 RandomTriangularMatrix (F, N, N, FFLAS::FflasUpper, FFLAS::FflasNonUnit, true, U, N, G); 487 // L <- U^T x R 488 FFLAS::fzero(F, N, N, L, N); 489 for (size_t k=0; k<R; ++k){ 490 size_t i = RRP[k]; 491 size_t j = CRP[k]; 492 FFLAS::fassign (F, N-i, U+i*(N+1), 1, L+j+i*N, N); 493 } 494 495 FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, N,N,N, F.one, L, N, U, N, F.zero, A, lda); 496 497 FFLAS::fflas_delete(L); 498 FFLAS::fflas_delete(U); 499 return A; 500 } 501 502 /** @brief Random Symmetric Matrix with prescribed rank and rank profile matrix 503 * Creates an \c n x \c n symmetric matrix with random entries and rank \c r. 504 * @param F field 505 * @param n order of \p A 506 * @param r rank of \p A 507 * @param A the matrix (preallocated to at least \c n x \c lda field elements) 508 * @param lda leading dimension of \p A 509 * @param RRP the R dimensional array with row positions of the rank profile matrix' pivots 510 * @param CRP the R dimensional array with column positions of the rank profile matrix' pivots 511 * @return \c A. 512 */ 513 template<class Field> 514 inline typename Field::Element_ptr RandomSymmetricMatrixWithRankandRPM(const Field & F,size_t M,size_t N,size_t R,typename Field::Element_ptr A,size_t lda,const size_t * RRP,const size_t * CRP)515 RandomSymmetricMatrixWithRankandRPM (const Field& F, size_t M, size_t N, size_t R, 516 typename Field::Element_ptr A, size_t lda, 517 const size_t * RRP, const size_t * CRP){ 518 typename Field::RandIter G(F); 519 return RandomSymmetricMatrixWithRankandRPM (F, N, R, A, lda, RRP, CRP, G); 520 } 521 522 /** @brief Random Matrix with prescribed rank, with random rank profile matrix 523 * Creates an \c m x \c n matrix with random entries, rank \c r and with a 524 * rank profile matrix chosen uniformly at random. 525 * @param F field 526 * @param m number of rows in \p A 527 * @param n number of cols in \p A 528 * @param r rank of the matrix to build 529 * @param A the matrix (preallocated to at least \c m x \c lda field elements) 530 * @param lda leading dimension of \p A 531 * @param G a random iterator 532 * @return \c A. 533 */ 534 template<class Field, class RandIter> 535 inline typename Field::Element_ptr RandomMatrixWithRankandRandomRPM(const Field & F,size_t M,size_t N,size_t R,typename Field::Element_ptr A,size_t lda,RandIter & G)536 RandomMatrixWithRankandRandomRPM (const Field& F, size_t M, size_t N, size_t R, 537 typename Field::Element_ptr A, size_t lda, RandIter& G){ 538 // generate the r pivots in the rank profile matrix E 539 size_t * pivot_r = FFLAS::fflas_new<size_t> (R); 540 size_t * pivot_c = FFLAS::fflas_new<size_t> (R); 541 RandomRankProfileMatrix (M, N, R, pivot_r, pivot_c); 542 RandomMatrixWithRankandRPM (F, M, N, R, A, lda, pivot_r, pivot_c, G); 543 FFLAS::fflas_delete(pivot_r); 544 FFLAS::fflas_delete(pivot_c); 545 return A; 546 } 547 548 /** @brief Random Matrix with prescribed rank, with random rank profile matrix 549 * Creates an \c m x \c n matrix with random entries, rank \c r and with a 550 * rank profile matrix chosen uniformly at random. 551 * @param F field 552 * @param m number of rows in \p A 553 * @param n number of cols in \p A 554 * @param r rank of the matrix to build 555 * @param A the matrix (preallocated to at least \c m x \c lda field elements) 556 * @param lda leading dimension of \p A 557 * @return \c A. 558 */ 559 template<class Field> 560 inline typename Field::Element_ptr RandomMatrixWithRankandRandomRPM(const Field & F,size_t M,size_t N,size_t R,typename Field::Element_ptr A,size_t lda)561 RandomMatrixWithRankandRandomRPM (const Field& F, size_t M, size_t N, size_t R, 562 typename Field::Element_ptr A, size_t lda){ 563 typename Field::RandIter G(F); 564 return RandomMatrixWithRankandRandomRPM (F, M, N, R, A, lda, G); 565 } 566 567 /** @brief Random Symmetric Matrix with prescribed rank, with random rank profile matrix 568 * Creates an \c n x \c n matrix with random entries, rank \c r and with a 569 * rank profile matrix chosen uniformly at random. 570 * @param F field 571 * @param n order of \p A 572 * @param r rank of \p A 573 * @param A the matrix (preallocated to at least \c n x \c lda field elements) 574 * @param lda leading dimension of \p A 575 * @param G a random iterator 576 * @return \c A. 577 */ 578 template<class Field, class RandIter> 579 inline typename Field::Element_ptr RandomSymmetricMatrixWithRankandRandomRPM(const Field & F,size_t N,size_t R,typename Field::Element_ptr A,size_t lda,RandIter & G)580 RandomSymmetricMatrixWithRankandRandomRPM (const Field& F, size_t N, size_t R, 581 typename Field::Element_ptr A, size_t lda, RandIter& G){ 582 // generate the r pivots in the rank profile matrix E 583 size_t * pivot_r = FFLAS::fflas_new<size_t> (R); 584 size_t * pivot_c = FFLAS::fflas_new<size_t> (R); 585 RandomSymmetricRankProfileMatrix (N, R, pivot_r, pivot_c); 586 RandomSymmetricMatrixWithRankandRPM (F, N, R, A, lda, pivot_r, pivot_c, G); 587 FFLAS::fflas_delete(pivot_r); 588 FFLAS::fflas_delete(pivot_c); 589 return A; 590 } 591 592 /** @brief Random Symmetric Matrix with prescribed rank, with random rank profile matrix 593 * Creates an \c n x \c n matrix with random entries, rank \c r and with a 594 * rank profile matrix chosen uniformly at random. 595 * @param F field 596 * @param n order of \p A 597 * @param r rank of \p A 598 * @param A the matrix (preallocated to at least \c n x \c lda field elements) 599 * @param lda leading dimension of \p A 600 * @return \c A. 601 */ 602 template<class Field> 603 inline typename Field::Element_ptr RandomSymmetricMatrixWithRankandRandomRPM(const Field & F,size_t N,size_t R,typename Field::Element_ptr A,size_t lda)604 RandomSymmetricMatrixWithRankandRandomRPM (const Field& F, size_t N, size_t R, 605 typename Field::Element_ptr A, size_t lda){ 606 typename Field::RandIter G(F); 607 return RandomSymmetricMatrixWithRankandRandomRPM (F, N, R, A, lda, G); 608 } 609 610 /** @brief Random Matrix with prescribed det. 611 * Creates a \c m x \c n matrix with random entries and rank \c r. 612 * @param F field 613 * @param d the prescribed value for the determinant of A 614 * @param n number of cols in \p A 615 * @param A the matrix to be generated (preallocated to at least \c n x \c lda field elements) 616 * @param lda leading dimension of \p A 617 * @return \c A. 618 */ 619 template<class Field> 620 inline typename Field::Element_ptr RandomMatrixWithDet(const Field & F,size_t n,const typename Field::Element d,typename Field::Element_ptr A,size_t lda)621 RandomMatrixWithDet(const Field & F, size_t n, const typename Field::Element d, 622 typename Field::Element_ptr A, size_t lda) { 623 typename Field::RandIter G(F); 624 return RandomMatrixWithDet (F, n, d, A, lda, G); 625 } 626 /** @brief Random Matrix with prescribed det. 627 * Creates a \c m x \c n matrix with random entries and rank \c r. 628 * @param F field 629 * @param d the prescribed value for the determinant of A 630 * @param n number of cols in \p A 631 * @param A the matrix to be generated (preallocated to at least \c n x \c lda field elements) 632 * @param lda leading dimension of \p A 633 * @return \c A. 634 */ 635 template<class Field, class RandIter> 636 inline typename Field::Element_ptr RandomMatrixWithDet(const Field & F,size_t n,const typename Field::Element d,typename Field::Element_ptr A,size_t lda,RandIter & G)637 RandomMatrixWithDet(const Field & F, size_t n, const typename Field::Element d, 638 typename Field::Element_ptr A, size_t lda, RandIter& G){ 639 FFLASFFPACK_check(n <= lda); 640 typedef typename Field::Element Element ; 641 Givaro::GeneralRingNonZeroRandIter<Field,RandIter> nzG (G); 642 643 size_t * P = FFLAS::fflas_new<size_t>(n); 644 size_t * Q = FFLAS::fflas_new<size_t>(n); 645 for (size_t i = 0 ; i < n ; ++i ) Q[i] = 0; 646 for (size_t i = 0 ; i < n ; ++i ) P[i] = 0; 647 648 Element * U = FFLAS::fflas_new<Element>(n*n); 649 Element * L = FFLAS::fflas_new<Element>(n*n); 650 651 /* Create a random P,Q */ 652 RandomPermutation (n, P); 653 RandomPermutation (n, Q); 654 655 /* det of P,Q */ 656 int d1 =1 ; 657 for (size_t i = 0 ; i < n ; ++i) 658 if (P[i] != i) 659 d1 = -d1; 660 for (size_t i = 0 ; i < n ; ++i) 661 if (Q[i] != i) 662 d1 = -d1; 663 664 /* Create L, lower det d */ 665 RandomTriangularMatrix (F, n, n, FFLAS::FflasLower, FFLAS::FflasNonUnit, true, L, n, G); 666 667 Element dd = F.one; 668 for (size_t i=0 ; i<n-1 ; ++i) 669 F.mulin(dd,L[i*n+i]); 670 671 F.div(dd,d,dd); 672 if (d1<0) F.negin(dd); 673 F.assign (L[n*n-1],dd); 674 675 /* Create U, upper unit*/ 676 RandomTriangularMatrix (F, n, n, FFLAS::FflasUpper, FFLAS::FflasUnit, true, U, n, G); 677 678 /* compute product */ 679 FFPACK::applyP (F, FFLAS::FflasRight, FFLAS::FflasNoTrans, 680 n,0,(int)n, U, n, P); 681 FFPACK::applyP (F, FFLAS::FflasLeft, FFLAS::FflasNoTrans, 682 n,0,(int)n, L, n, Q); 683 FFLAS::fgemm (F, FFLAS::FflasNoTrans, FFLAS::FflasNoTrans, 684 n,n,n, 1.0, L, n, U, n, 0.0, A, lda); 685 // @todo compute LU with ftrtr 686 687 FFLAS::fflas_delete( P); 688 FFLAS::fflas_delete( L); 689 FFLAS::fflas_delete( U); 690 FFLAS::fflas_delete( Q); 691 692 return A; 693 } 694 } // FFPACK 695 #endif 696 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 697 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 698