1 /* initRandom.c */
2
3 #include "../../Drand.h"
4 #include "../SubMtx.h"
5
6 /*--------------------------------------------------------------------*/
7 /*
8 ------------------------------------------------
9 purpose -- to initialize a matrix object with
10 random entries and possibly random structure.
11
12 created -- 98feb16, cca
13 ------------------------------------------------
14 */
15 void
SubMtx_initRandom(SubMtx * mtx,int type,int mode,int rowid,int colid,int nrow,int ncol,int nent,int seed)16 SubMtx_initRandom (
17 SubMtx *mtx,
18 int type,
19 int mode,
20 int rowid,
21 int colid,
22 int nrow,
23 int ncol,
24 int nent,
25 int seed
26 ) {
27 Drand *drand ;
28 /*
29 ---------------
30 check the input
31 ---------------
32 */
33 if ( mtx == NULL || type < 1 || type > 2 || mode < 0 || mode > 9
34 || nrow <= 0 || ncol <= 0 ) {
35 fprintf(stderr, "\n fatal error in SubMtx_initRandom()"
36 "\n bad input\n") ;
37 exit(-1) ;
38 }
39 SubMtx_clearData(mtx) ;
40 /*
41 --------------------------------------
42 initialize the random number generator
43 --------------------------------------
44 */
45 drand = Drand_new() ;
46 if ( seed > 0 ) {
47 Drand_setSeed(drand, seed) ;
48 }
49 Drand_setUniform(drand, 0, 1) ;
50 /*
51 ---------------------------
52 switch over the matrix mode
53 ---------------------------
54 */
55 switch ( mode ) {
56 case SUBMTX_DENSE_ROWS :
57 case SUBMTX_DENSE_COLUMNS : {
58 double *entries ;
59 int inc1, inc2 ;
60 /*
61 ----------------------------------------
62 this case is simple, fill entries vector
63 ----------------------------------------
64 */
65 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nrow*ncol) ;
66 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
67 if ( SUBMTX_IS_REAL(mtx) ) {
68 Drand_fillDvector(drand, nrow*ncol, entries) ;
69 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
70 Drand_fillDvector(drand, 2*nrow*ncol, entries) ;
71 }
72 } break ;
73 case SUBMTX_SPARSE_ROWS :
74 case SUBMTX_SPARSE_COLUMNS :
75 case SUBMTX_SPARSE_TRIPLES : {
76 double *entries ;
77 int ii ;
78 int *colids, *indices, *ivec, *ivec1, *ivec2, *rowids, *sizes ;
79 /*
80 ----------------------------------------------------------
81 (1) fill ivec[nrow*ncol] with 0, 1, ..., nrow*ncol-1
82 (2) shuffle ivec[nrow*ncol]
83 (3) set rowids[nent] to be the row id
84 of the entries of ivec[0:nent-1]
85 (4) set colids[nent] to be the column id
86 of the entries of ivec[0:nent-1]
87 (5) for sparse rows and columns, set sizes[] and indices[]
88 ----------------------------------------------------------
89 */
90 ivec = IVinit(nrow*ncol, -1) ;
91 IVramp(nrow*ncol, ivec, 0, 1) ;
92 IVshuffle(nrow*ncol, ivec, seed + 3) ;
93 if ( nent > nrow*ncol ) {
94 nent = nrow*ncol ;
95 }
96 rowids = IVinit(nent, -1) ;
97 colids = IVinit(nent, -1) ;
98 for ( ii = 0 ; ii < nent ; ii++ ) {
99 rowids[ii] = ivec[ii] % nrow ;
100 colids[ii] = ivec[ii] / nrow ;
101 }
102 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
103 switch ( mode ) {
104 case SUBMTX_SPARSE_ROWS :
105 IV2qsortUp(nent, rowids, colids) ;
106 SubMtx_sparseRowsInfo(mtx, &nrow, &nent,
107 &sizes, &indices, &entries) ;
108 IVzero(nrow, sizes) ;
109 for ( ii = 0 ; ii < nent ; ii++ ) {
110 sizes[rowids[ii]]++ ;
111 }
112 IVcopy(nent, indices, colids) ;
113 break ;
114 case SUBMTX_SPARSE_COLUMNS :
115 IV2qsortUp(nent, colids, rowids) ;
116 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
117 &sizes, &indices, &entries) ;
118 IVzero(ncol, sizes) ;
119 for ( ii = 0 ; ii < nent ; ii++ ) {
120 sizes[colids[ii]]++ ;
121 }
122 IVcopy(nent, indices, rowids) ;
123 break ;
124 case SUBMTX_SPARSE_TRIPLES :
125 IV2qsortUp(nent, colids, rowids) ;
126 SubMtx_sparseTriplesInfo(mtx, &nent, &ivec1, &ivec2, &entries) ;
127 IVcopy(nent, ivec1, rowids) ;
128 IVcopy(nent, ivec2, colids) ;
129 break ;
130 }
131 if ( SUBMTX_IS_REAL(mtx) ) {
132 Drand_fillDvector(drand, nent, entries) ;
133 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
134 Drand_fillDvector(drand, 2*nent, entries) ;
135 }
136 IVfree(ivec) ;
137 IVfree(rowids) ;
138 IVfree(colids) ;
139 } break ;
140 case SUBMTX_DENSE_SUBROWS : {
141 double *entries ;
142 int irow, size, *firstlocs, *ivec1, *ivec2, *sizes ;
143
144 ivec1 = IVinit(nrow, -1) ;
145 ivec2 = IVinit(nrow, 0) ;
146 Drand_setUniform(drand, 0, ncol) ;
147 Drand_fillIvector(drand, nrow, ivec1) ;
148 Drand_fillIvector(drand, nrow, ivec2) ;
149 for ( irow = nent = 0 ; irow < nrow ; irow++ ) {
150 if ( (size = ivec2[irow] - ivec1[irow] + 1) <= 0 ) {
151 ivec1[irow] = -1 ;
152 ivec2[irow] = 0 ;
153 } else {
154 nent += size ;
155 ivec2[irow] = size ;
156 }
157 }
158 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
159 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
160 &firstlocs, &sizes, &entries) ;
161 IVcopy(nrow, firstlocs, ivec1) ;
162 IVcopy(nrow, sizes, ivec2) ;
163 if ( SUBMTX_IS_REAL(mtx) ) {
164 Drand_fillDvector(drand, nent, entries) ;
165 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
166 Drand_fillDvector(drand, 2*nent, entries) ;
167 }
168 IVfree(ivec1) ;
169 IVfree(ivec2) ;
170 } break ;
171 case SUBMTX_DENSE_SUBCOLUMNS : {
172 double *entries ;
173 int jcol, size, *firstlocs, *ivec1, *ivec2, *sizes ;
174
175 ivec1 = IVinit(ncol, -1) ;
176 ivec2 = IVinit(ncol, -1) ;
177 Drand_setUniform(drand, 0, nrow) ;
178 Drand_fillIvector(drand, ncol, ivec1) ;
179 Drand_fillIvector(drand, ncol, ivec2) ;
180 for ( jcol = nent = 0 ; jcol < ncol ; jcol++ ) {
181 if ( (size = ivec2[jcol] - ivec1[jcol] + 1) <= 0 ) {
182 ivec1[jcol] = -1 ;
183 ivec2[jcol] = 0 ;
184 } else {
185 nent += size ;
186 ivec2[jcol] = size ;
187 }
188 }
189 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
190 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
191 &firstlocs, &sizes, &entries) ;
192 IVcopy(ncol, firstlocs, ivec1) ;
193 IVcopy(ncol, sizes, ivec2) ;
194 if ( SUBMTX_IS_REAL(mtx) ) {
195 Drand_fillDvector(drand, nent, entries) ;
196 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
197 Drand_fillDvector(drand, 2*nent, entries) ;
198 }
199 IVfree(ivec1) ;
200 IVfree(ivec2) ;
201 } break ;
202 case SUBMTX_DIAGONAL : {
203 double *entries ;
204
205 ncol = nrow ;
206 nent = nrow ;
207 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
208 SubMtx_diagonalInfo(mtx, &ncol, &entries) ;
209 if ( SUBMTX_IS_REAL(mtx) ) {
210 Drand_fillDvector(drand, nent, entries) ;
211 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
212 Drand_fillDvector(drand, 2*nent, entries) ;
213 }
214 } break ;
215 case SUBMTX_BLOCK_DIAGONAL_SYM : {
216 double *entries ;
217 int ipivot, irow, m ;
218 int *pivots, *pivotsizes ;
219
220 ncol = nrow ;
221 pivotsizes = IVinit(nrow, -1) ;
222 Drand_setUniform(drand, 1,3) ;
223 Drand_fillIvector(drand, ncol, pivotsizes) ;
224 for ( ipivot = irow = nent = 0 ; irow < nrow ; ipivot++ ) {
225 m = pivotsizes[ipivot] ;
226 if ( irow + m > nrow ) {
227 m = pivotsizes[ipivot] = nrow - irow ;
228 }
229 nent += (m*(m+1))/2 ;
230 irow += m ;
231 }
232 IVzero(nrow - ipivot, pivotsizes + ipivot) ;
233 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
234 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent, &pivots, &entries) ;
235 IVcopy(nrow, pivots, pivotsizes) ;
236 if ( SUBMTX_IS_REAL(mtx) ) {
237 Drand_fillDvector(drand, nent, entries) ;
238 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
239 Drand_fillDvector(drand, 2*nent, entries) ;
240 }
241 IVfree(pivotsizes) ;
242 } break ;
243 case SUBMTX_BLOCK_DIAGONAL_HERM : {
244 double *entries ;
245 int ipivot, irow, kk, m ;
246 int *pivots, *pivotsizes ;
247
248 ncol = nrow ;
249 pivotsizes = IVinit(nrow, -1) ;
250 Drand_setUniform(drand, 1,3) ;
251 Drand_fillIvector(drand, ncol, pivotsizes) ;
252 for ( ipivot = irow = nent = 0 ; irow < nrow ; ipivot++ ) {
253 m = pivotsizes[ipivot] ;
254 if ( irow + m > nrow ) {
255 m = pivotsizes[ipivot] = nrow - irow ;
256 }
257 nent += (m*(m+1))/2 ;
258 irow += m ;
259 }
260 IVzero(nrow - ipivot, pivotsizes + ipivot) ;
261 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
262 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent, &pivots, &entries) ;
263 IVcopy(nrow, pivots, pivotsizes) ;
264 if ( SUBMTX_IS_REAL(mtx) ) {
265 Drand_fillDvector(drand, nent, entries) ;
266 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
267 Drand_fillDvector(drand, 2*nent, entries) ;
268 }
269 if ( SUBMTX_IS_COMPLEX(mtx) ) {
270 for ( ipivot = irow = kk = 0 ; irow < nrow ; ipivot++ ) {
271 m = pivotsizes[ipivot] ;
272 if ( m == 1 ) {
273 entries[2*kk+1] = 0.0 ;
274 kk++ ;
275 } else {
276 entries[2*kk+1] = 0.0 ;
277 entries[2*kk+5] = 0.0 ;
278 kk += 3 ;
279 }
280 irow += m ;
281 }
282 }
283 IVfree(pivotsizes) ;
284 } break ;
285 default :
286 fprintf(stderr, "\n\n %% type %d not yet supported", type) ;
287 exit(0) ;
288 break ;
289 }
290 Drand_free(drand) ;
291
292 return ; }
293
294 /*--------------------------------------------------------------------*/
295 /*
296 ------------------------------------------------
297 purpose -- to initialize a matrix object with
298 random entries in the upper triangle
299
300 strict = 1 --> strict upper triangle
301
302 created -- 98feb16, cca
303 ------------------------------------------------
304 */
305 void
SubMtx_initRandomUpperTriangle(SubMtx * mtx,int type,int mode,int rowid,int colid,int nrow,int ncol,int nent,int seed,int strict)306 SubMtx_initRandomUpperTriangle (
307 SubMtx *mtx,
308 int type,
309 int mode,
310 int rowid,
311 int colid,
312 int nrow,
313 int ncol,
314 int nent,
315 int seed,
316 int strict
317 ) {
318 Drand *drand ;
319 /*
320 ---------------
321 check the input
322 ---------------
323 */
324 if ( mtx == NULL || type < 1 || type > 2 || mode < 0 || mode > 9
325 || nrow <= 0 || ncol <= 0 ) {
326 fprintf(stderr, "\n fatal error in SubMtx_initRandomUpperTriangle()"
327 "\n bad input\n") ;
328 exit(-1) ;
329 }
330 SubMtx_clearData(mtx) ;
331 /*
332 --------------------------------------
333 initialize the random number generator
334 --------------------------------------
335 */
336 drand = Drand_new() ;
337 if ( seed > 0 ) {
338 Drand_setSeed(drand, seed) ;
339 }
340 Drand_setUniform(drand, 0, 1) ;
341 /*
342 ---------------------------
343 switch over the matrix type
344 ---------------------------
345 */
346 switch ( mode ) {
347 case SUBMTX_DENSE_ROWS :
348 case SUBMTX_DENSE_COLUMNS : {
349 double *entries ;
350 int ij, inc1, inc2, irow, jcol ;
351 /*
352 -------------------------------------------------------------------
353 (1) initialize the matrix and fill all entries with random numbers.
354 (2) put zeros in (strict) lower triangle
355 -------------------------------------------------------------------
356 */
357 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nrow*ncol) ;
358 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
359 if ( SUBMTX_IS_REAL(mtx) ) {
360 Drand_fillDvector(drand, nent, entries) ;
361 if ( strict == 1 ) {
362 for ( irow = 0 ; irow < nrow ; irow++ ) {
363 for ( jcol = 0 ; jcol <= irow ; jcol++ ) {
364 ij = irow*inc1 + jcol*inc2 ;
365 entries[ij] = 0.0 ;
366 }
367 }
368 } else {
369 for ( irow = 0 ; irow < nrow ; irow++ ) {
370 for ( jcol = 0 ; jcol < irow ; jcol++ ) {
371 ij = irow*inc1 + jcol*inc2 ;
372 entries[ij] = 0.0 ;
373 }
374 }
375 }
376 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
377 Drand_fillDvector(drand, 2*nent, entries) ;
378 if ( strict == 1 ) {
379 for ( irow = 0 ; irow < nrow ; irow++ ) {
380 for ( jcol = 0 ; jcol <= irow ; jcol++ ) {
381 ij = irow*inc1 + jcol*inc2 ;
382 entries[2*ij] = entries[2*ij+1] = 0.0 ;
383 }
384 }
385 } else {
386 for ( irow = 0 ; irow < nrow ; irow++ ) {
387 for ( jcol = 0 ; jcol < irow ; jcol++ ) {
388 ij = irow*inc1 + jcol*inc2 ;
389 entries[2*ij] = entries[2*ij+1] = 0.0 ;
390 }
391 }
392 }
393 }
394 } break ;
395 case SUBMTX_SPARSE_ROWS :
396 case SUBMTX_SPARSE_COLUMNS : {
397 double *entries ;
398 int count, ii, irow, jcol ;
399 int *colids, *indices, *ivec, *rowids, *sizes ;
400 /*
401 ------------------------------------------
402 fill ivec[*] with upper triangular entries
403 ------------------------------------------
404 */
405 ivec = IVinit(nrow*ncol, -1) ;
406 if ( strict == 1 ) {
407 for ( irow = count = 0 ; irow < nrow ; irow++ ) {
408 for ( jcol = irow + 1 ; jcol < ncol ; jcol++ ) {
409 ivec[count++] = irow + jcol*nrow ;
410 }
411 }
412 } else {
413 for ( irow = count = 0 ; irow < nrow ; irow++ ) {
414 for ( jcol = irow ; jcol < ncol ; jcol++ ) {
415 ivec[count++] = irow + jcol*nrow ;
416 }
417 }
418 }
419 IVshuffle(count, ivec, seed) ;
420 if ( nent > count ) {
421 nent = count ;
422 }
423 rowids = IVinit(nent, -1) ;
424 colids = IVinit(nent, -1) ;
425 for ( ii = 0 ; ii < nent ; ii++ ) {
426 rowids[ii] = ivec[ii] % nrow ;
427 colids[ii] = ivec[ii] / nrow ;
428 }
429 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
430 if ( mode == SUBMTX_SPARSE_ROWS ) {
431 SubMtx_sparseRowsInfo(mtx, &nrow, &nent,
432 &sizes, &indices, &entries) ;
433 IVzero(nrow, sizes) ;
434 IV2qsortUp(nent, rowids, colids) ;
435 for ( ii = 0 ; ii < nent ; ii++ ) {
436 sizes[rowids[ii]]++ ;
437 }
438 IVcopy(nent, indices, colids) ;
439 } else {
440 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
441 &sizes, &indices, &entries) ;
442 IV2qsortUp(nent, colids, rowids) ;
443 IVzero(ncol, sizes) ;
444 for ( ii = 0 ; ii < nent ; ii++ ) {
445 sizes[colids[ii]]++ ;
446 }
447 IVcopy(nent, indices, rowids) ;
448 }
449 if ( SUBMTX_IS_REAL(mtx) ) {
450 Drand_fillDvector(drand, nent, entries) ;
451 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
452 Drand_fillDvector(drand, 2*nent, entries) ;
453 }
454 IVfree(ivec) ;
455 IVfree(rowids) ;
456 IVfree(colids) ;
457 } break ;
458 case SUBMTX_DENSE_SUBROWS : {
459 double *entries ;
460 int irow, size ;
461 int *firstlocs, *ivec1, *ivec2, *sizes ;
462 /*
463 ------------------------------------------------------------
464 fill ivec1[] and ivec2[] with start and end column locations
465 ------------------------------------------------------------
466 */
467 ivec1 = IVinit(nrow, -1) ;
468 ivec2 = IVinit(nrow, -1) ;
469 Drand_setUniform(drand, 0, ncol) ;
470 Drand_fillIvector(drand, nrow, ivec1) ;
471 Drand_fillIvector(drand, nrow, ivec2) ;
472 /*
473 -------------------------------------------------
474 modify ivec1[] to be on or just past the diagonal
475 -------------------------------------------------
476 */
477 if ( strict == 1 ) {
478 for ( irow = 0 ; irow < nrow ; irow++ ) {
479 if ( ivec1[irow] <= irow ) {
480 ivec1[irow] = irow + 1 ;
481 }
482 }
483 } else {
484 for ( irow = 0 ; irow < nrow ; irow++ ) {
485 if ( ivec1[irow] < irow ) {
486 ivec1[irow] = irow ;
487 }
488 }
489 }
490 /*
491 -----------------------------------------------------------
492 set ivec2[] to be the sizes and count the number of entries
493 -----------------------------------------------------------
494 */
495 for ( irow = nent = 0 ; irow < nrow ; irow++ ) {
496 if ( (size = ivec2[irow] - ivec1[irow] + 1) <= 0 ) {
497 ivec1[irow] = -1 ;
498 ivec2[irow] = 0 ;
499 } else {
500 nent += size ;
501 ivec2[irow] = size ;
502 }
503 }
504 /*
505 -------------------------------------------------
506 initialize, set first location, sizes and entries
507 -------------------------------------------------
508 */
509 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
510 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
511 &firstlocs, &sizes, &entries) ;
512 IVcopy(nrow, firstlocs, ivec1) ;
513 IVcopy(nrow, sizes, ivec2) ;
514 if ( SUBMTX_IS_REAL(mtx) ) {
515 Drand_fillDvector(drand, nent, entries) ;
516 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
517 Drand_fillDvector(drand, 2*nent, entries) ;
518 }
519 IVfree(ivec1) ;
520 IVfree(ivec2) ;
521 } break ;
522 case SUBMTX_DENSE_SUBCOLUMNS : {
523 double *entries ;
524 int jcol, size ;
525 int *firstlocs, *ivec1, *ivec2, *sizes ;
526 /*
527 ---------------------------------------------------------
528 fill ivec1[] and ivec2[] with start and end row locations
529 ---------------------------------------------------------
530 */
531 ivec1 = IVinit(ncol, -1) ;
532 ivec2 = IVinit(ncol, -1) ;
533 Drand_setUniform(drand, 0, nrow) ;
534 Drand_fillIvector(drand, ncol, ivec1) ;
535 Drand_fillIvector(drand, ncol, ivec2) ;
536 /*
537 -------------------------------------------------
538 modify ivec2[] to be on or just past the diagonal
539 -------------------------------------------------
540 */
541 if ( strict == 1 ) {
542 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
543 if ( ivec2[jcol] >= jcol ) {
544 ivec2[jcol] = jcol - 1 ;
545 }
546 }
547 } else {
548 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
549 if ( ivec2[jcol] > jcol ) {
550 ivec2[jcol] = jcol ;
551 }
552 }
553 }
554 /*
555 -----------------------------------------------------------
556 set ivec2[] to be the sizes and count the number of entries
557 -----------------------------------------------------------
558 */
559 for ( jcol = nent = 0 ; jcol < ncol ; jcol++ ) {
560 if ( (size = ivec2[jcol] - ivec1[jcol] + 1) <= 0 ) {
561 ivec1[jcol] = -1 ;
562 ivec2[jcol] = 0 ;
563 } else {
564 nent += size ;
565 ivec2[jcol] = size ;
566 }
567 }
568 /*
569 -------------------------------------------------
570 initialize, set first location, sizes and entries
571 -------------------------------------------------
572 */
573 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
574 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
575 &firstlocs, &sizes, &entries) ;
576 IVcopy(nrow, firstlocs, ivec1) ;
577 IVcopy(nrow, sizes, ivec2) ;
578 if ( SUBMTX_IS_REAL(mtx) ) {
579 Drand_fillDvector(drand, nent, entries) ;
580 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
581 Drand_fillDvector(drand, 2*nent, entries) ;
582 }
583 IVfree(ivec1) ;
584 IVfree(ivec2) ;
585 } break ;
586 default :
587 fprintf(stderr, "\n\n %% type %d not yet supported", type) ;
588 exit(0) ;
589 break ;
590 }
591 /*
592 --------------------------------
593 free the random number generator
594 --------------------------------
595 */
596 Drand_free(drand) ;
597
598 return ; }
599
600 /*--------------------------------------------------------------------*/
601 /*
602 ------------------------------------------------
603 purpose -- to initialize a matrix object with
604 random entries in the lower triangle
605
606 strict = 1 --> strict lower triangle
607
608 created -- 98feb16, cca
609 ------------------------------------------------
610 */
611 void
SubMtx_initRandomLowerTriangle(SubMtx * mtx,int type,int mode,int rowid,int colid,int nrow,int ncol,int nent,int seed,int strict)612 SubMtx_initRandomLowerTriangle (
613 SubMtx *mtx,
614 int type,
615 int mode,
616 int rowid,
617 int colid,
618 int nrow,
619 int ncol,
620 int nent,
621 int seed,
622 int strict
623 ) {
624 Drand *drand ;
625 int *colind, *rowind ;
626 /*
627 ---------------
628 check the input
629 ---------------
630 */
631 if ( mtx == NULL || type < 0 || type > 9 || nrow <= 0 || ncol <= 0 ) {
632 fprintf(stderr, "\n fatal error in SubMtx_initRandomLowerTriangle()"
633 "\n bad input\n") ;
634 exit(-1) ;
635 }
636 SubMtx_clearData(mtx) ;
637 /*
638 --------------------------------------
639 initialize the random number generator
640 --------------------------------------
641 */
642 drand = Drand_new() ;
643 if ( seed > 0 ) {
644 Drand_setSeed(drand, seed) ;
645 }
646 Drand_setUniform(drand, 0, 1) ;
647 /*
648 ---------------------------
649 switch over the matrix type
650 ---------------------------
651 */
652 switch ( mode ) {
653 case SUBMTX_DENSE_ROWS :
654 case SUBMTX_DENSE_COLUMNS : {
655 double *entries ;
656 int ij, inc1, inc2, irow, jcol ;
657 /*
658 -------------------------------------------------------------------
659 (1) initialize the matrix and fill all entries with random numbers.
660 (2) put zeros in (strict) upper triangle
661 -------------------------------------------------------------------
662 */
663 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nrow*ncol) ;
664 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
665 if ( SUBMTX_IS_REAL(mtx) ) {
666 Drand_fillDvector(drand, nent, entries) ;
667 if ( strict == 1 ) {
668 for ( irow = 0 ; irow < nrow ; irow++ ) {
669 for ( jcol = irow ; jcol < ncol ; jcol++ ) {
670 ij = irow*inc1 + jcol*inc2 ;
671 entries[ij] = 0.0 ;
672 }
673 }
674 } else {
675 for ( irow = 0 ; irow < nrow ; irow++ ) {
676 for ( jcol = irow + 1 ; jcol < ncol ; jcol++ ) {
677 ij = irow*inc1 + jcol*inc2 ;
678 entries[ij] = 0.0 ;
679 }
680 }
681 }
682 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
683 Drand_fillDvector(drand, 2*nent, entries) ;
684 if ( strict == 1 ) {
685 for ( irow = 0 ; irow < nrow ; irow++ ) {
686 for ( jcol = irow ; jcol < ncol ; jcol++ ) {
687 ij = irow*inc1 + jcol*inc2 ;
688 entries[2*ij] = entries[2*ij+1] = 0.0 ;
689 }
690 }
691 } else {
692 for ( irow = 0 ; irow < nrow ; irow++ ) {
693 for ( jcol = irow + 1 ; jcol < ncol ; jcol++ ) {
694 ij = irow*inc1 + jcol*inc2 ;
695 entries[2*ij] = entries[2*ij+1] = 0.0 ;
696 }
697 }
698 }
699 }
700 } break ;
701 case SUBMTX_SPARSE_ROWS :
702 case SUBMTX_SPARSE_COLUMNS : {
703 double *entries ;
704 int count, ii, irow, jcol ;
705 int *colids, *indices, *ivec, *rowids, *sizes ;
706
707 ivec = IVinit(nrow*ncol, -1) ;
708 if ( strict == 1 ) {
709 for ( irow = count = 0 ; irow < nrow ; irow++ ) {
710 for ( jcol = 0 ; jcol < irow ; jcol++ ) {
711 ivec[count++] = irow + jcol*nrow ;
712 }
713 }
714 } else {
715 for ( irow = count = 0 ; irow < nrow ; irow++ ) {
716 for ( jcol = 0 ; jcol <= irow ; jcol++ ) {
717 ivec[count++] = irow + jcol*nrow ;
718 }
719 }
720 }
721 IVshuffle(count, ivec, seed) ;
722 if ( nent > count ) {
723 nent = count ;
724 }
725 rowids = IVinit(nent, -1) ;
726 colids = IVinit(nent, -1) ;
727 for ( ii = 0 ; ii < nent ; ii++ ) {
728 rowids[ii] = ivec[ii] % nrow ;
729 colids[ii] = ivec[ii] / nrow ;
730 }
731 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
732 if ( mode == SUBMTX_SPARSE_ROWS ) {
733 SubMtx_sparseRowsInfo(mtx, &nrow, &nent,
734 &sizes, &indices, &entries) ;
735 IVzero(nrow, sizes) ;
736 IV2qsortUp(nent, rowids, colids) ;
737 for ( ii = 0 ; ii < nent ; ii++ ) {
738 sizes[rowids[ii]]++ ;
739 }
740 IVcopy(nent, indices, colids) ;
741 } else {
742 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
743 &sizes, &indices, &entries) ;
744 IV2qsortUp(nent, colids, rowids) ;
745 IVzero(ncol, sizes) ;
746 for ( ii = 0 ; ii < nent ; ii++ ) {
747 sizes[colids[ii]]++ ;
748 }
749 IVcopy(nent, indices, rowids) ;
750 }
751 if ( SUBMTX_IS_REAL(mtx) ) {
752 Drand_fillDvector(drand, nent, entries) ;
753 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
754 Drand_fillDvector(drand, 2*nent, entries) ;
755 }
756 IVfree(ivec) ;
757 IVfree(rowids) ;
758 IVfree(colids) ;
759 } break ;
760 case SUBMTX_DENSE_SUBROWS : {
761 double *entries ;
762 int irow, size ;
763 int *firstlocs, *ivec1, *ivec2, *sizes ;
764 /*
765 ------------------------------------------------------------
766 fill ivec1[] and ivec2[] with start and end column locations
767 ------------------------------------------------------------
768 */
769 ivec1 = IVinit(nrow, -1) ;
770 ivec2 = IVinit(nrow, -1) ;
771 Drand_setUniform(drand, 0, ncol) ;
772 Drand_fillIvector(drand, nrow, ivec1) ;
773 Drand_fillIvector(drand, nrow, ivec2) ;
774 /*
775 -------------------------------------------------
776 modify ivec2[] to be on or just past the diagonal
777 -------------------------------------------------
778 */
779 if ( strict == 1 ) {
780 for ( irow = 0 ; irow < nrow ; irow++ ) {
781 if ( ivec2[irow] >= irow ) {
782 ivec2[irow] = irow - 1 ;
783 }
784 }
785 } else {
786 for ( irow = 0 ; irow < nrow ; irow++ ) {
787 if ( ivec1[irow] > irow ) {
788 ivec1[irow] = irow ;
789 }
790 }
791 }
792 /*
793 -----------------------------------------------------------
794 set ivec2[] to be the sizes and count the number of entries
795 -----------------------------------------------------------
796 */
797 for ( irow = nent = 0 ; irow < nrow ; irow++ ) {
798 if ( (size = ivec2[irow] - ivec1[irow] + 1) <= 0 ) {
799 ivec1[irow] = -1 ;
800 ivec2[irow] = 0 ;
801 } else {
802 nent += size ;
803 ivec2[irow] = size ;
804 }
805 }
806 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
807 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
808 &firstlocs, &sizes, &entries) ;
809 IVcopy(nrow, firstlocs, ivec1) ;
810 IVcopy(nrow, sizes, ivec2) ;
811 if ( SUBMTX_IS_REAL(mtx) ) {
812 Drand_fillDvector(drand, nent, entries) ;
813 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
814 Drand_fillDvector(drand, 2*nent, entries) ;
815 }
816 IVfree(ivec1) ;
817 IVfree(ivec2) ;
818 } break ;
819 case SUBMTX_DENSE_SUBCOLUMNS : {
820 double *entries ;
821 int jcol, size ;
822 int *firstlocs, *ivec1, *ivec2, *sizes ;
823 /*
824 ---------------------------------------------------------
825 fill ivec1[] and ivec2[] with start and end row locations
826 ---------------------------------------------------------
827 */
828 ivec1 = IVinit(ncol, -1) ;
829 ivec2 = IVinit(ncol, -1) ;
830 Drand_setUniform(drand, 0, nrow) ;
831 Drand_fillIvector(drand, ncol, ivec1) ;
832 Drand_fillIvector(drand, ncol, ivec2) ;
833 /*
834 -------------------------------------------------
835 modify ivec1[] to be on or just past the diagonal
836 -------------------------------------------------
837 */
838 if ( strict == 1 ) {
839 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
840 if ( ivec1[jcol] <= jcol ) {
841 ivec1[jcol] = jcol + 1 ;
842 }
843 }
844 } else {
845 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
846 if ( ivec1[jcol] < jcol ) {
847 ivec1[jcol] = jcol ;
848 }
849 }
850 }
851 /*
852 -----------------------------------------------------------
853 set ivec2[] to be the sizes and count the number of entries
854 -----------------------------------------------------------
855 */
856 for ( jcol = nent = 0 ; jcol < ncol ; jcol++ ) {
857 if ( (size = ivec2[jcol] - ivec1[jcol] + 1) <= 0 ) {
858 ivec1[jcol] = -1 ;
859 ivec2[jcol] = 0 ;
860 } else {
861 nent += size ;
862 }
863 }
864 SubMtx_init(mtx, type, mode, rowid, colid, nrow, ncol, nent) ;
865 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
866 &firstlocs, &sizes, &entries) ;
867 IVcopy(ncol, firstlocs, ivec1) ;
868 IVcopy(ncol, sizes, ivec2) ;
869 if ( SUBMTX_IS_REAL(mtx) ) {
870 Drand_fillDvector(drand, nent, entries) ;
871 } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
872 Drand_fillDvector(drand, 2*nent, entries) ;
873 }
874 IVfree(ivec1) ;
875 IVfree(ivec2) ;
876 } break ;
877 default :
878 fprintf(stderr, "\n\n %% type %d not yet supported", type) ;
879 exit(0) ;
880 break ;
881 }
882 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
883 IVramp(nrow, rowind, 0, 1) ;
884 SubMtx_columnIndices(mtx, &ncol, &colind) ;
885 IVramp(ncol, colind, 0, 1) ;
886 /*
887 --------------------------------
888 free the random number generator
889 --------------------------------
890 */
891 Drand_free(drand) ;
892
893 return ; }
894
895 /*--------------------------------------------------------------------*/
896