SCM

SCM Repository

[matrix] Annotation of /pkg/Matrix/src/CHOLMOD/Core/cholmod_transpose.c
ViewVC logotype

Annotation of /pkg/Matrix/src/CHOLMOD/Core/cholmod_transpose.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2798 - (view) (download) (as text)

1 : bates 912 /* ========================================================================== */
2 :     /* === Core/cholmod_transpose =============================================== */
3 :     /* ========================================================================== */
4 :    
5 :     /* -----------------------------------------------------------------------------
6 : bates 1876 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 : bates 1515 * Univ. of Florida. Author: Timothy A. Davis
8 : bates 912 * The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
9 :     * Lesser General Public License. See lesser.txt for a text of the license.
10 :     * CHOLMOD is also available under other licenses; contact authors for details.
11 :     * -------------------------------------------------------------------------- */
12 :    
13 :     /* Core utility routines for the cholmod_sparse object to
14 :     * compute the transpose or permuted transpose of a matrix:
15 :     *
16 :     * Primary routines:
17 :     * -----------------
18 :     * cholmod_transpose transpose sparse matrix
19 :     * cholmod_ptranspose transpose and permute sparse matrix
20 :     * cholmod_sort sort row indices in each column of sparse matrix
21 :     *
22 :     * Secondary routines:
23 :     * -------------------
24 :     * cholmod_transpose_unsym transpose unsymmetric sparse matrix
25 :     * cholmod_transpose_sym transpose symmetric sparse matrix
26 :     *
27 :     * All xtypes (pattern, real, complex, and zomplex) are supported.
28 :     *
29 :     * ---------------------------------------
30 :     * Unsymmetric case: A->stype is zero.
31 :     * ---------------------------------------
32 :     *
33 :     * Computes F = A', F = A (:,f)' or F = A (p,f)', except that the indexing by
34 :     * f does not work the same as the MATLAB notation (see below). A->stype
35 :     * is zero, which denotes that both the upper and lower triangular parts of
36 :     * A are present (and used). A may in fact be symmetric in pattern and/or
37 :     * value; A->stype just denotes which part of A are stored. A may be
38 :     * rectangular.
39 :     *
40 :     * p is a permutation of 0:m-1, and f is a subset of 0:n-1, where A is m-by-n.
41 :     * There can be no duplicate entries in p or f.
42 :     *
43 :     * The set f is held in fset and fsize.
44 : bates 1515 * fset = NULL means ":" in MATLAB. fsize is ignored.
45 :     * fset != NULL means f = fset [0..fsize-1].
46 : bates 912 * fset != NULL and fsize = 0 means f is the empty set.
47 :     *
48 :     * Columns not in the set f are considered to be zero. That is,
49 :     * if A is 5-by-10 then F = A (:,[3 4])' is not 2-by-5, but 10-by-5, and rows
50 :     * 3 and 4 of F are equal to columns 3 and 4 of A (the other rows of F are
51 :     * zero). More precisely, in MATLAB notation:
52 :     *
53 :     * [m n] = size (A) ;
54 :     * F = A ;
55 :     * notf = ones (1,n) ;
56 :     * notf (f) = 0 ;
57 :     * F (:, find (notf)) = 0
58 :     * F = F'
59 :     *
60 :     * If you want the MATLAB equivalent F=A(p,f) operation, use cholmod_submatrix
61 :     * instead (which does not compute the transpose).
62 :     *
63 :     * F->nzmax must be large enough to hold the matrix F. It is not modified.
64 :     * If F->nz is present then F->nz [j] = # of entries in column j of F.
65 :     *
66 :     * A can be sorted or unsorted, with packed or unpacked columns.
67 :     *
68 :     * If f is present and not sorted in ascending order, then F is unsorted
69 :     * (that is, it may contain columns whose row indices do not appear in
70 :     * ascending order). Otherwise, F is sorted (the row indices in each
71 :     * column of F appear in strictly ascending order).
72 :     *
73 :     * F is returned in packed or unpacked form, depending on F->packed on input.
74 :     * If F->packed is false, then F is returned in unpacked form (F->nz must be
75 :     * present). Each row i of F is large enough to hold all the entries in row i
76 :     * of A, even if f is provided. That is, F->i and
77 :     * F->x [F->p [i] .. F->p [i] + F->nz [i] - 1] contain all entries in A (i,f),
78 :     * but F->p [i+1] - F->p [i] is equal to the number of nonzeros in A (i,:),
79 :     * not just A (i,f).
80 :     *
81 :     * The cholmod_transpose_unsym routine is the only operation in CHOLMOD that
82 :     * can produce an unpacked matrix.
83 :     *
84 :     * ---------------------------------------
85 :     * Symmetric case: A->stype is nonzero.
86 :     * ---------------------------------------
87 :     *
88 :     * Computes F = A' or F = A(p,p)', the transpose or permuted transpose, where
89 :     * A->stype is nonzero.
90 :     *
91 :     * If A->stype > 0, then A is a symmetric matrix where just the upper part
92 :     * of the matrix is stored. Entries in the lower triangular part may be
93 :     * present, but are ignored. A must be square. If F=A', then F is returned
94 :     * sorted; otherwise F is unsorted for the F=A(p,p)' case.
95 :     *
96 :     * There can be no duplicate entries in p.
97 :     * The fset and fsize parameters are not used.
98 :     *
99 :     * Three kinds of transposes are available, depending on the "values" parameter:
100 :     * 0: do not transpose the numerical values; create a CHOLMOD_PATTERN matrix
101 :     * 1: array transpose
102 :     * 2: complex conjugate transpose (same as 2 if input is real or pattern)
103 :     *
104 :     * -----------------------------------------------------------------------------
105 :     *
106 :     * For cholmod_transpose_unsym and cholmod_transpose_sym, the output matrix
107 :     * F must already be pre-allocated by the caller, with the correct dimensions.
108 :     * If F is not valid or has the wrong dimensions, it is not modified.
109 :     * Otherwise, if F is too small, the transpose is not computed; the contents
110 :     * of F->p contain the column pointers of the resulting matrix, where
111 :     * F->p [F->ncol] > F->nzmax. In this case, the remaining contents of F are
112 :     * not modified. F can still be properly free'd with cholmod_free_sparse.
113 :     */
114 :    
115 : bates 1515 #include "cholmod_internal.h"
116 : bates 912 #include "cholmod_core.h"
117 :    
118 :    
119 :     /* ========================================================================== */
120 :     /* === TEMPLATE ============================================================= */
121 :     /* ========================================================================== */
122 :    
123 :     #define PATTERN
124 :     #include "t_cholmod_transpose.c"
125 :     #define REAL
126 :     #include "t_cholmod_transpose.c"
127 :     #define COMPLEX
128 :     #include "t_cholmod_transpose.c"
129 :     #define COMPLEX
130 :     #define NCONJUGATE
131 :     #include "t_cholmod_transpose.c"
132 :     #define ZOMPLEX
133 :     #include "t_cholmod_transpose.c"
134 :     #define ZOMPLEX
135 :     #define NCONJUGATE
136 :     #include "t_cholmod_transpose.c"
137 :    
138 :    
139 :     /* ========================================================================== */
140 :     /* === cholmod_transpose_unsym ============================================== */
141 :     /* ========================================================================== */
142 :    
143 :     /* Compute F = A', A (:,f)', or A (p,f)', where A is unsymmetric and F is
144 :     * already allocated. See cholmod_transpose for a simpler routine.
145 :     *
146 :     * workspace:
147 :     * Iwork (MAX (nrow,ncol)) if fset is present
148 :     * Iwork (nrow) if fset is NULL
149 :     *
150 :     * The xtype of A and F must match, unless values is zero or F->xtype is
151 :     * CHOLMOD_PATTERN (in which case only the pattern of A is transpose into F).
152 :     */
153 :    
154 :     int CHOLMOD(transpose_unsym)
155 :     (
156 :     /* ---- input ---- */
157 :     cholmod_sparse *A, /* matrix to transpose */
158 :     int values, /* 2: complex conj. transpose, 1: array transpose,
159 :     0: do not transpose the numerical values */
160 :     Int *Perm, /* size nrow, if present (can be NULL) */
161 :     Int *fset, /* subset of 0:(A->ncol)-1 */
162 :     size_t fsize, /* size of fset */
163 :     /* ---- output --- */
164 :     cholmod_sparse *F, /* F = A', A(:,f)', or A(p,f)' */
165 :     /* --------------- */
166 :     cholmod_common *Common
167 :     )
168 :     {
169 :     Int *Fp, *Fnz, *Ap, *Ai, *Anz, *Wi ;
170 : bates 1515 Int nrow, ncol, permute, use_fset, Apacked, Fpacked, p, pend,
171 : bates 912 i, j, k, Fsorted, nf, jj, jlast ;
172 : bates 1515 size_t s ;
173 :     int ok = TRUE ;
174 : bates 912
175 :     /* ---------------------------------------------------------------------- */
176 :     /* check inputs */
177 :     /* ---------------------------------------------------------------------- */
178 :    
179 :     RETURN_IF_NULL_COMMON (FALSE) ;
180 :     RETURN_IF_NULL (A, FALSE) ;
181 :     RETURN_IF_NULL (F, FALSE) ;
182 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
183 :     RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
184 :     if (A->nrow != F->ncol || A->ncol != F->nrow)
185 :     {
186 :     ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
187 :     return (FALSE) ;
188 :     }
189 :     Common->status = CHOLMOD_OK ;
190 :    
191 :     /* ---------------------------------------------------------------------- */
192 :     /* get inputs */
193 :     /* ---------------------------------------------------------------------- */
194 :    
195 :     nf = fsize ;
196 :     use_fset = (fset != NULL) ;
197 :     nrow = A->nrow ;
198 :     ncol = A->ncol ;
199 :    
200 :     Ap = A->p ; /* size A->ncol+1, column pointers of A */
201 :     Ai = A->i ; /* size nz = Ap [A->ncol], row indices of A */
202 :     Anz = A->nz ;
203 :     Apacked = A->packed ;
204 :     ASSERT (IMPLIES (!Apacked, Anz != NULL)) ;
205 :    
206 :     permute = (Perm != NULL) ;
207 :    
208 :     Fp = F->p ; /* size A->nrow+1, row pointers of F */
209 :     Fnz = F->nz ;
210 :     Fpacked = F->packed ;
211 :     ASSERT (IMPLIES (!Fpacked, Fnz != NULL)) ;
212 :    
213 :     nf = (use_fset) ? nf : ncol ;
214 :    
215 :     /* ---------------------------------------------------------------------- */
216 :     /* allocate workspace */
217 :     /* ---------------------------------------------------------------------- */
218 :    
219 : bates 1515 /* s = nrow + ((fset != NULL) ? ncol : 0) */
220 :     s = CHOLMOD(add_size_t) (nrow, ((fset != NULL) ? ncol : 0), &ok) ;
221 :     if (!ok)
222 :     {
223 :     ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
224 :     return (FALSE) ;
225 :     }
226 :    
227 :     CHOLMOD(allocate_work) (0, s, 0, Common) ;
228 : bates 912 if (Common->status < CHOLMOD_OK)
229 :     {
230 :     return (FALSE) ; /* out of memory */
231 :     }
232 :    
233 :     Wi = Common->Iwork ; /* size nrow (i/l/l) */
234 :    
235 :     /* ---------------------------------------------------------------------- */
236 :     /* check Perm and fset */
237 :     /* ---------------------------------------------------------------------- */
238 :    
239 :     if (permute)
240 :     {
241 :     for (i = 0 ; i < nrow ; i++)
242 :     {
243 :     Wi [i] = 1 ;
244 :     }
245 :     for (k = 0 ; k < nrow ; k++)
246 :     {
247 :     i = Perm [k] ;
248 :     if (i < 0 || i > nrow || Wi [i] == 0)
249 :     {
250 :     ERROR (CHOLMOD_INVALID, "invalid permutation") ;
251 :     return (FALSE) ;
252 :     }
253 :     Wi [i] = 0 ;
254 :     }
255 :     }
256 :    
257 :     if (use_fset)
258 :     {
259 :     for (j = 0 ; j < ncol ; j++)
260 :     {
261 :     Wi [j] = 1 ;
262 :     }
263 :     for (k = 0 ; k < nf ; k++)
264 :     {
265 :     j = fset [k] ;
266 :     if (j < 0 || j > ncol || Wi [j] == 0)
267 :     {
268 :     ERROR (CHOLMOD_INVALID, "invalid fset") ;
269 :     return (FALSE) ;
270 :     }
271 :     Wi [j] = 0 ;
272 :     }
273 :     }
274 :    
275 :     /* Perm and fset are now valid */
276 :     ASSERT (CHOLMOD(dump_perm) (Perm, nrow, nrow, "Perm", Common)) ;
277 :     ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;
278 :    
279 :     /* ---------------------------------------------------------------------- */
280 :     /* count the entries in each row of A or A(:,f) */
281 :     /* ---------------------------------------------------------------------- */
282 :    
283 :     for (i = 0 ; i < nrow ; i++)
284 :     {
285 :     Wi [i] = 0 ;
286 :     }
287 :    
288 :     jlast = EMPTY ;
289 :     Fsorted = TRUE ;
290 :    
291 :     if (use_fset)
292 :     {
293 :     /* count entries in each row of A(:,f) */
294 :     for (jj = 0 ; jj < nf ; jj++)
295 :     {
296 :     j = fset [jj] ;
297 :     if (j <= jlast)
298 :     {
299 :     Fsorted = FALSE ;
300 :     }
301 :     p = Ap [j] ;
302 :     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
303 :     for ( ; p < pend ; p++)
304 :     {
305 :     Wi [Ai [p]]++ ;
306 :     }
307 :     jlast = j ;
308 :     }
309 :    
310 :     /* save the nz counts if F is unpacked, and recount all of A */
311 :     if (!Fpacked)
312 :     {
313 :     if (permute)
314 :     {
315 :     for (i = 0 ; i < nrow ; i++)
316 :     {
317 :     Fnz [i] = Wi [Perm [i]] ;
318 :     }
319 :     }
320 :     else
321 :     {
322 :     for (i = 0 ; i < nrow ; i++)
323 :     {
324 :     Fnz [i] = Wi [i] ;
325 :     }
326 :     }
327 :     for (i = 0 ; i < nrow ; i++)
328 :     {
329 :     Wi [i] = 0 ;
330 :     }
331 :    
332 :     /* count entries in each row of A */
333 :     for (j = 0 ; j < ncol ; j++)
334 :     {
335 :     p = Ap [j] ;
336 :     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
337 :     for ( ; p < pend ; p++)
338 :     {
339 :     Wi [Ai [p]]++ ;
340 :     }
341 :     }
342 :     }
343 :    
344 :     }
345 :     else
346 :     {
347 :    
348 :     /* count entries in each row of A */
349 :     for (j = 0 ; j < ncol ; j++)
350 :     {
351 :     p = Ap [j] ;
352 :     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
353 :     for ( ; p < pend ; p++)
354 :     {
355 :     Wi [Ai [p]]++ ;
356 :     }
357 :     }
358 :    
359 :     /* save the nz counts if F is unpacked */
360 :     if (!Fpacked)
361 :     {
362 :     if (permute)
363 :     {
364 :     for (i = 0 ; i < nrow ; i++)
365 :     {
366 :     Fnz [i] = Wi [Perm [i]] ;
367 :     }
368 :     }
369 :     else
370 :     {
371 :     for (i = 0 ; i < nrow ; i++)
372 :     {
373 :     Fnz [i] = Wi [i] ;
374 :     }
375 :     }
376 :     }
377 :     }
378 :    
379 :     /* ---------------------------------------------------------------------- */
380 :     /* compute the row pointers */
381 :     /* ---------------------------------------------------------------------- */
382 :    
383 :     p = 0 ;
384 :     if (permute)
385 :     {
386 :     for (i = 0 ; i < nrow ; i++)
387 :     {
388 :     Fp [i] = p ;
389 :     p += Wi [Perm [i]] ;
390 :     }
391 :     for (i = 0 ; i < nrow ; i++)
392 :     {
393 :     Wi [Perm [i]] = Fp [i] ;
394 :     }
395 :     }
396 :     else
397 :     {
398 :     for (i = 0 ; i < nrow ; i++)
399 :     {
400 :     Fp [i] = p ;
401 :     p += Wi [i] ;
402 :     }
403 :     for (i = 0 ; i < nrow ; i++)
404 :     {
405 :     Wi [i] = Fp [i] ;
406 :     }
407 :     }
408 :     Fp [nrow] = p ;
409 :    
410 :     if (p > (Int) (F->nzmax))
411 :     {
412 :     ERROR (CHOLMOD_INVALID, "F is too small") ;
413 :     return (FALSE) ;
414 :     }
415 :    
416 :     /* ---------------------------------------------------------------------- */
417 :     /* transpose matrix, using template routine */
418 :     /* ---------------------------------------------------------------------- */
419 :    
420 : bates 1515 ok = FALSE ;
421 : bates 912 if (values == 0 || F->xtype == CHOLMOD_PATTERN)
422 :     {
423 :     ok = p_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
424 :     }
425 :     else if (F->xtype == CHOLMOD_REAL)
426 :     {
427 :     ok = r_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
428 :     }
429 :     else if (F->xtype == CHOLMOD_COMPLEX)
430 :     {
431 :     if (values == 1)
432 :     {
433 :     /* array transpose */
434 :     ok = ct_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
435 :     }
436 :     else
437 :     {
438 :     /* complex conjugate transpose */
439 :     ok = c_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
440 :     }
441 :     }
442 :     else if (F->xtype == CHOLMOD_ZOMPLEX)
443 :     {
444 :     if (values == 1)
445 :     {
446 :     /* array transpose */
447 :     ok = zt_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
448 :     }
449 :     else
450 :     {
451 :     /* complex conjugate transpose */
452 :     ok = z_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
453 :     }
454 :     }
455 :    
456 :     /* ---------------------------------------------------------------------- */
457 :     /* finalize result F */
458 :     /* ---------------------------------------------------------------------- */
459 :    
460 :     if (ok)
461 :     {
462 :     F->sorted = Fsorted ;
463 :     }
464 :     ASSERT (CHOLMOD(dump_sparse) (F, "output F unsym", Common) >= 0) ;
465 :     return (ok) ;
466 :     }
467 :    
468 :    
469 :     /* ========================================================================== */
470 :     /* === cholmod_transpose_sym ================================================ */
471 :     /* ========================================================================== */
472 :    
473 :     /* Compute F = A' or A (p,p)', where A is symmetric and F is already allocated.
474 :     * See cholmod_transpose for a simpler routine.
475 :     *
476 :     * workspace: Iwork (nrow) if Perm NULL, Iwork (2*nrow) if Perm non-NULL.
477 :     */
478 :    
479 :     int CHOLMOD(transpose_sym)
480 :     (
481 :     /* ---- input ---- */
482 :     cholmod_sparse *A, /* matrix to transpose */
483 :     int values, /* 2: complex conj. transpose, 1: array transpose,
484 :     0: do not transpose the numerical values */
485 :     Int *Perm, /* size nrow, if present (can be NULL) */
486 :     /* ---- output --- */
487 :     cholmod_sparse *F, /* F = A' or A(p,p)' */
488 :     /* --------------- */
489 :     cholmod_common *Common
490 :     )
491 :     {
492 :     Int *Ap, *Anz, *Ai, *Fp, *Wi, *Pinv, *Iwork ;
493 :     Int p, pend, packed, upper, permute, jold, n, i, j, k, iold ;
494 : bates 1515 size_t s ;
495 :     int ok = TRUE ;
496 : bates 912
497 :     /* ---------------------------------------------------------------------- */
498 :     /* check inputs */
499 :     /* ---------------------------------------------------------------------- */
500 :    
501 :     RETURN_IF_NULL_COMMON (FALSE) ;
502 :     RETURN_IF_NULL (A, FALSE) ;
503 :     RETURN_IF_NULL (F, FALSE) ;
504 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
505 :     RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
506 :     if (A->nrow != A->ncol || A->stype == 0)
507 :     {
508 :     /* this routine handles square symmetric matrices only */
509 :     ERROR (CHOLMOD_INVALID, "matrix must be symmetric") ;
510 :     return (FALSE) ;
511 :     }
512 :     if (A->nrow != F->ncol || A->ncol != F->nrow)
513 :     {
514 :     ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
515 :     return (FALSE) ;
516 :     }
517 :     Common->status = CHOLMOD_OK ;
518 :    
519 :     /* ---------------------------------------------------------------------- */
520 :     /* get inputs */
521 :     /* ---------------------------------------------------------------------- */
522 :    
523 :     permute = (Perm != NULL) ;
524 :     n = A->nrow ;
525 :     Ap = A->p ; /* size A->ncol+1, column pointers of A */
526 :     Ai = A->i ; /* size nz = Ap [A->ncol], row indices of A */
527 :     Anz = A->nz ;
528 :     packed = A->packed ;
529 :     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
530 :     upper = (A->stype > 0) ;
531 :    
532 :     Fp = F->p ; /* size A->nrow+1, row pointers of F */
533 :    
534 :     /* ---------------------------------------------------------------------- */
535 :     /* allocate workspace */
536 :     /* ---------------------------------------------------------------------- */
537 :    
538 : bates 1515 /* s = (Perm != NULL) ? 2*n : n */
539 :     s = CHOLMOD(add_size_t) (n, ((Perm != NULL) ? n : 0), &ok) ;
540 :     if (!ok)
541 :     {
542 :     ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
543 :     return (FALSE) ;
544 :     }
545 :    
546 :     CHOLMOD(allocate_work) (0, s, 0, Common) ;
547 : bates 912 if (Common->status < CHOLMOD_OK)
548 :     {
549 :     return (FALSE) ; /* out of memory */
550 :     }
551 :    
552 :     /* ---------------------------------------------------------------------- */
553 :     /* get workspace */
554 :     /* ---------------------------------------------------------------------- */
555 :    
556 :     Iwork = Common->Iwork ;
557 :     Wi = Iwork ; /* size n (i/l/l) */
558 :     Pinv = Iwork + n ; /* size n (i/i/l) , unused if Perm NULL */
559 :    
560 :     /* ---------------------------------------------------------------------- */
561 :     /* check Perm and construct inverse permutation */
562 :     /* ---------------------------------------------------------------------- */
563 :    
564 :     if (permute)
565 :     {
566 :     for (i = 0 ; i < n ; i++)
567 :     {
568 :     Pinv [i] = EMPTY ;
569 :     }
570 :     for (k = 0 ; k < n ; k++)
571 :     {
572 :     i = Perm [k] ;
573 :     if (i < 0 || i > n || Pinv [i] != EMPTY)
574 :     {
575 :     ERROR (CHOLMOD_INVALID, "invalid permutation") ;
576 :     return (FALSE) ;
577 :     }
578 :     Pinv [i] = k ;
579 :     }
580 :     }
581 :    
582 :     /* Perm is now valid */
583 :     ASSERT (CHOLMOD(dump_perm) (Perm, n, n, "Perm", Common)) ;
584 :    
585 :     /* ---------------------------------------------------------------------- */
586 :     /* count the entries in each row of F */
587 :     /* ---------------------------------------------------------------------- */
588 :    
589 :     for (i = 0 ; i < n ; i++)
590 :     {
591 :     Wi [i] = 0 ;
592 :     }
593 :    
594 :     if (packed)
595 :     {
596 :     if (permute)
597 :     {
598 :     if (upper)
599 :     {
600 :     /* packed, permuted, upper */
601 :     for (j = 0 ; j < n ; j++)
602 :     {
603 :     jold = Perm [j] ;
604 :     pend = Ap [jold+1] ;
605 :     for (p = Ap [jold] ; p < pend ; p++)
606 :     {
607 :     iold = Ai [p] ;
608 :     if (iold <= jold)
609 :     {
610 :     i = Pinv [iold] ;
611 :     Wi [MIN (i, j)]++ ;
612 :     }
613 :     }
614 :     }
615 :     }
616 :     else
617 :     {
618 :     /* packed, permuted, lower */
619 :     for (j = 0 ; j < n ; j++)
620 :     {
621 :     jold = Perm [j] ;
622 :     pend = Ap [jold+1] ;
623 :     for (p = Ap [jold] ; p < pend ; p++)
624 :     {
625 :     iold = Ai [p] ;
626 :     if (iold >= jold)
627 :     {
628 :     i = Pinv [iold] ;
629 :     Wi [MAX (i, j)]++ ;
630 :     }
631 :     }
632 :     }
633 :     }
634 :     }
635 :     else
636 :     {
637 :     if (upper)
638 :     {
639 :     /* packed, unpermuted, upper */
640 :     for (j = 0 ; j < n ; j++)
641 :     {
642 :     pend = Ap [j+1] ;
643 :     for (p = Ap [j] ; p < pend ; p++)
644 :     {
645 :     i = Ai [p] ;
646 :     if (i <= j)
647 :     {
648 :     Wi [i]++ ;
649 :     }
650 :     }
651 :     }
652 :     }
653 :     else
654 :     {
655 :     /* packed, unpermuted, lower */
656 :     for (j = 0 ; j < n ; j++)
657 :     {
658 :     pend = Ap [j+1] ;
659 :     for (p = Ap [j] ; p < pend ; p++)
660 :     {
661 :     i = Ai [p] ;
662 :     if (i >= j)
663 :     {
664 :     Wi [i]++ ;
665 :     }
666 :     }
667 :     }
668 :     }
669 :     }
670 :     }
671 :     else
672 :     {
673 :     if (permute)
674 :     {
675 :     if (upper)
676 :     {
677 :     /* unpacked, permuted, upper */
678 :     for (j = 0 ; j < n ; j++)
679 :     {
680 :     jold = Perm [j] ;
681 :     p = Ap [jold] ;
682 :     pend = p + Anz [jold] ;
683 :     for ( ; p < pend ; p++)
684 :     {
685 :     iold = Ai [p] ;
686 :     if (iold <= jold)
687 :     {
688 :     i = Pinv [iold] ;
689 :     Wi [MIN (i, j)]++ ;
690 :     }
691 :     }
692 :     }
693 :     }
694 :     else
695 :     {
696 :     /* unpacked, permuted, lower */
697 :     for (j = 0 ; j < n ; j++)
698 :     {
699 :     jold = Perm [j] ;
700 :     p = Ap [jold] ;
701 :     pend = p + Anz [jold] ;
702 :     for ( ; p < pend ; p++)
703 :     {
704 :     iold = Ai [p] ;
705 :     if (iold >= jold)
706 :     {
707 :     i = Pinv [iold] ;
708 :     Wi [MAX (i, j)]++ ;
709 :     }
710 :     }
711 :     }
712 :     }
713 :     }
714 :     else
715 :     {
716 :     if (upper)
717 :     {
718 :     /* unpacked, unpermuted, upper */
719 :     for (j = 0 ; j < n ; j++)
720 :     {
721 :     p = Ap [j] ;
722 :     pend = p + Anz [j] ;
723 :     for ( ; p < pend ; p++)
724 :     {
725 :     i = Ai [p] ;
726 :     if (i <= j)
727 :     {
728 :     Wi [i]++ ;
729 :     }
730 :     }
731 :     }
732 :     }
733 :     else
734 :     {
735 :     /* unpacked, unpermuted, lower */
736 :     for (j = 0 ; j < n ; j++)
737 :     {
738 :     p = Ap [j] ;
739 :     pend = p + Anz [j] ;
740 :     for ( ; p < pend ; p++)
741 :     {
742 :     i = Ai [p] ;
743 :     if (i >= j)
744 :     {
745 :     Wi [i]++ ;
746 :     }
747 :     }
748 :     }
749 :     }
750 :     }
751 :     }
752 :    
753 :     /* ---------------------------------------------------------------------- */
754 :     /* compute the row pointers */
755 :     /* ---------------------------------------------------------------------- */
756 :    
757 :     p = 0 ;
758 :     for (i = 0 ; i < n ; i++)
759 :     {
760 :     Fp [i] = p ;
761 :     p += Wi [i] ;
762 :     }
763 :     Fp [n] = p ;
764 :     for (i = 0 ; i < n ; i++)
765 :     {
766 :     Wi [i] = Fp [i] ;
767 :     }
768 :    
769 :     if (p > (Int) (F->nzmax))
770 :     {
771 :     ERROR (CHOLMOD_INVALID, "F is too small") ;
772 :     return (FALSE) ;
773 :     }
774 :    
775 :     /* ---------------------------------------------------------------------- */
776 :     /* transpose matrix, using template routine */
777 :     /* ---------------------------------------------------------------------- */
778 :    
779 : bates 1515 ok = FALSE ;
780 : bates 912 if (values == 0 || F->xtype == CHOLMOD_PATTERN)
781 :     {
782 :     PRINT2 (("\n:::: p_transpose_sym Perm %p\n", Perm)) ;
783 :     ok = p_cholmod_transpose_sym (A, Perm, F, Common) ;
784 :     }
785 :     else if (F->xtype == CHOLMOD_REAL)
786 :     {
787 :     PRINT2 (("\n:::: r_transpose_sym Perm %p\n", Perm)) ;
788 :     ok = r_cholmod_transpose_sym (A, Perm, F, Common) ;
789 :     }
790 :     else if (F->xtype == CHOLMOD_COMPLEX)
791 :     {
792 :     if (values == 1)
793 :     {
794 :     /* array transpose */
795 :     PRINT2 (("\n:::: ct_transpose_sym Perm %p\n", Perm)) ;
796 :     ok = ct_cholmod_transpose_sym (A, Perm, F, Common) ;
797 :     }
798 :     else
799 :     {
800 :     /* complex conjugate transpose */
801 :     PRINT2 (("\n:::: c_transpose_sym Perm %p\n", Perm)) ;
802 :     ok = c_cholmod_transpose_sym (A, Perm, F, Common) ;
803 :     }
804 :     }
805 :     else if (F->xtype == CHOLMOD_ZOMPLEX)
806 :     {
807 :     if (values == 1)
808 :     {
809 :     /* array transpose */
810 :     PRINT2 (("\n:::: zt_transpose_sym Perm %p\n", Perm)) ;
811 :     ok = zt_cholmod_transpose_sym (A, Perm, F, Common) ;
812 :     }
813 :     else
814 :     {
815 :     /* complex conjugate transpose */
816 :     PRINT2 (("\n:::: z_transpose_sym Perm %p\n", Perm)) ;
817 :     ok = z_cholmod_transpose_sym (A, Perm, F, Common) ;
818 :     }
819 :     }
820 :    
821 :     /* ---------------------------------------------------------------------- */
822 :     /* finalize result F */
823 :     /* ---------------------------------------------------------------------- */
824 :    
825 :     /* F is sorted if there is no permutation vector */
826 :     if (ok)
827 :     {
828 :     F->sorted = !permute ;
829 :     F->packed = TRUE ;
830 :     F->stype = - SIGN (A->stype) ; /* flip the stype */
831 :     ASSERT (CHOLMOD(dump_sparse) (F, "output F sym", Common) >= 0) ;
832 :     }
833 :     return (ok) ;
834 :     }
835 :    
836 :    
837 :     /* ========================================================================== */
838 :     /* === cholmod_transpose ==================================================== */
839 :     /* ========================================================================== */
840 :    
841 :     /* Returns A'. See also cholmod_ptranspose below. */
842 :    
843 :     cholmod_sparse *CHOLMOD(transpose)
844 :     (
845 :     /* ---- input ---- */
846 :     cholmod_sparse *A, /* matrix to transpose */
847 :     int values, /* 2: complex conj. transpose, 1: array transpose,
848 :     0: do not transpose the numerical values
849 :     (returns its result as CHOLMOD_PATTERN) */
850 :     /* --------------- */
851 :     cholmod_common *Common
852 :     )
853 :     {
854 :     return (CHOLMOD(ptranspose) (A, values, NULL, NULL, 0, Common)) ;
855 :     }
856 :    
857 :    
858 :     /* ========================================================================== */
859 :     /* === cholmod_ptranspose =================================================== */
860 :     /* ========================================================================== */
861 :    
862 :     /* Return A' or A(p,p)' if A is symmetric. Return A', A(:,f)', or A(p,f)' if
863 :     * A is unsymmetric.
864 :     *
865 :     * workspace:
866 :     * Iwork (MAX (nrow,ncol)) if unsymmetric and fset is non-NULL
867 :     * Iwork (nrow) if unsymmetric and fset is NULL
868 :     * Iwork (2*nrow) if symmetric and Perm is non-NULL.
869 :     * Iwork (nrow) if symmetric and Perm is NULL.
870 :     *
871 :     * A simple worst-case upper bound on the workspace is nrow+ncol.
872 :     */
873 :    
874 :     cholmod_sparse *CHOLMOD(ptranspose)
875 :     (
876 :     /* ---- input ---- */
877 :     cholmod_sparse *A, /* matrix to transpose */
878 :     int values, /* 2: complex conj. transpose, 1: array transpose,
879 :     0: do not transpose the numerical values */
880 :     Int *Perm, /* if non-NULL, F = A(p,f) or A(p,p) */
881 :     Int *fset, /* subset of 0:(A->ncol)-1 */
882 :     size_t fsize, /* size of fset */
883 :     /* --------------- */
884 :     cholmod_common *Common
885 :     )
886 :     {
887 :     Int *Ap, *Anz ;
888 :     cholmod_sparse *F ;
889 : bates 1515 Int nrow, ncol, use_fset, j, jj, fnz, packed, stype, nf, xtype ;
890 :     size_t ineed ;
891 :     int ok = TRUE ;
892 : bates 912
893 :     nf = fsize ;
894 :    
895 :     /* ---------------------------------------------------------------------- */
896 :     /* check inputs */
897 :     /* ---------------------------------------------------------------------- */
898 :    
899 :     RETURN_IF_NULL_COMMON (NULL) ;
900 :     RETURN_IF_NULL (A, FALSE) ;
901 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
902 :     stype = A->stype ;
903 :     Common->status = CHOLMOD_OK ;
904 :    
905 :     /* ---------------------------------------------------------------------- */
906 :     /* allocate workspace */
907 :     /* ---------------------------------------------------------------------- */
908 :    
909 :     nrow = A->nrow ;
910 :     ncol = A->ncol ;
911 :    
912 :     if (stype != 0)
913 :     {
914 :     use_fset = FALSE ;
915 :     if (Perm != NULL)
916 :     {
917 : bates 1515 ineed = CHOLMOD(mult_size_t) (A->nrow, 2, &ok) ;
918 : bates 912 }
919 :     else
920 :     {
921 : bates 1515 ineed = A->nrow ;
922 : bates 912 }
923 :     }
924 :     else
925 :     {
926 :     use_fset = (fset != NULL) ;
927 :     if (use_fset)
928 :     {
929 : bates 1515 ineed = MAX (A->nrow, A->ncol) ;
930 : bates 912 }
931 :     else
932 :     {
933 : bates 1515 ineed = A->nrow ;
934 : bates 912 }
935 :     }
936 :    
937 : bates 1515 if (!ok)
938 :     {
939 :     ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
940 :     return (NULL) ;
941 :     }
942 :    
943 : bates 912 CHOLMOD(allocate_work) (0, ineed, 0, Common) ;
944 :     if (Common->status < CHOLMOD_OK)
945 :     {
946 :     return (NULL) ; /* out of memory */
947 :     }
948 :    
949 :     /* ---------------------------------------------------------------------- */
950 :     /* get inputs */
951 :     /* ---------------------------------------------------------------------- */
952 :    
953 :     Ap = A->p ;
954 :     Anz = A->nz ;
955 :     packed = A->packed ;
956 :     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
957 :     xtype = values ? A->xtype : CHOLMOD_PATTERN ;
958 :    
959 :     /* ---------------------------------------------------------------------- */
960 :     /* allocate F */
961 :     /* ---------------------------------------------------------------------- */
962 :    
963 :     /* determine # of nonzeros in F */
964 :     if (stype != 0)
965 :     {
966 :     /* F=A' or F=A(p,p)', fset is ignored */
967 :     fnz = CHOLMOD(nnz) (A, Common) ;
968 :     }
969 :     else
970 :     {
971 :     nf = (use_fset) ? nf : ncol ;
972 :     if (use_fset)
973 :     {
974 :     fnz = 0 ;
975 :     /* F=A(:,f)' or F=A(p,f)' */
976 :     for (jj = 0 ; jj < nf ; jj++)
977 :     {
978 :     /* The fset is not yet checked; it will be thoroughly checked
979 :     * in cholmod_transpose_unsym. For now, just make sure we don't
980 :     * access Ap and Anz out of bounds. */
981 :     j = fset [jj] ;
982 :     if (j >= 0 && j < ncol)
983 :     {
984 :     fnz += packed ? (Ap [j+1] - Ap [j]) : MAX (0, Anz [j]) ;
985 :     }
986 :     }
987 :     }
988 :     else
989 :     {
990 :     /* F=A' or F=A(p,:)' */
991 :     fnz = CHOLMOD(nnz) (A, Common) ;
992 :     }
993 :     }
994 :    
995 :     /* F is ncol-by-nrow, fnz nonzeros, sorted unless f is present and unsorted,
996 :     * packed, of opposite stype as A, and with/without numerical values */
997 :     F = CHOLMOD(allocate_sparse) (ncol, nrow, fnz, TRUE, TRUE, -SIGN(stype),
998 :     xtype, Common) ;
999 :     if (Common->status < CHOLMOD_OK)
1000 :     {
1001 :     return (NULL) ; /* out of memory */
1002 :     }
1003 :    
1004 :     /* ---------------------------------------------------------------------- */
1005 :     /* transpose and optionally permute the matrix A */
1006 :     /* ---------------------------------------------------------------------- */
1007 :    
1008 :     if (stype != 0)
1009 :     {
1010 :     /* F = A (p,p)', using upper or lower triangular part of A only */
1011 :     ok = CHOLMOD(transpose_sym) (A, values, Perm, F, Common) ;
1012 :     }
1013 :     else
1014 :     {
1015 :     /* F = A (p,f)' */
1016 :     ok = CHOLMOD(transpose_unsym) (A, values, Perm, fset, nf, F, Common) ;
1017 :     }
1018 :    
1019 :     /* ---------------------------------------------------------------------- */
1020 :     /* return the matrix F, or NULL if an error occured */
1021 :     /* ---------------------------------------------------------------------- */
1022 :    
1023 :     if (!ok)
1024 :     {
1025 :     CHOLMOD(free_sparse) (&F, Common) ;
1026 :     }
1027 :     return (F) ;
1028 :     }
1029 :    
1030 :    
1031 :     /* ========================================================================== */
1032 :     /* === cholmod_sort ========================================================= */
1033 :     /* ========================================================================== */
1034 :    
1035 :     /* Sort the columns of A, in place. Returns A in packed form, even if it
1036 :     * starts as unpacked. Removes entries in the ignored part of a symmetric
1037 :     * matrix.
1038 :     *
1039 :     * workspace: Iwork (max (nrow,ncol)). Allocates additional workspace for a
1040 :     * temporary copy of A'.
1041 :     */
1042 :    
1043 :     int CHOLMOD(sort)
1044 :     (
1045 :     /* ---- in/out --- */
1046 :     cholmod_sparse *A, /* matrix to sort */
1047 :     /* --------------- */
1048 :     cholmod_common *Common
1049 :     )
1050 :     {
1051 :     Int *Ap ;
1052 :     cholmod_sparse *F ;
1053 :     Int anz, ncol, nrow, stype ;
1054 :    
1055 :     /* ---------------------------------------------------------------------- */
1056 :     /* check inputs */
1057 :     /* ---------------------------------------------------------------------- */
1058 :    
1059 :     RETURN_IF_NULL_COMMON (FALSE) ;
1060 :     RETURN_IF_NULL (A, FALSE) ;
1061 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
1062 :     Common->status = CHOLMOD_OK ;
1063 :     nrow = A->nrow ;
1064 :     if (nrow <= 1)
1065 :     {
1066 :     /* a 1-by-n sparse matrix must be sorted */
1067 :     A->sorted = TRUE ;
1068 :     return (TRUE) ;
1069 :     }
1070 :    
1071 :     /* ---------------------------------------------------------------------- */
1072 :     /* allocate workspace */
1073 :     /* ---------------------------------------------------------------------- */
1074 :    
1075 :     ncol = A->ncol ;
1076 :     CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
1077 :     if (Common->status < CHOLMOD_OK)
1078 :     {
1079 :     return (FALSE) ; /* out of memory */
1080 :     }
1081 :    
1082 :     /* ---------------------------------------------------------------------- */
1083 :     /* get inputs */
1084 :     /* ---------------------------------------------------------------------- */
1085 :    
1086 :     anz = CHOLMOD(nnz) (A, Common) ;
1087 :     stype = A->stype ;
1088 :    
1089 :     /* ---------------------------------------------------------------------- */
1090 :     /* sort the columns of the matrix */
1091 :     /* ---------------------------------------------------------------------- */
1092 :    
1093 :     /* allocate workspace for transpose: ncol-by-nrow, same # of nonzeros as A,
1094 :     * sorted, packed, same stype as A, and of the same numeric type as A. */
1095 :     F = CHOLMOD(allocate_sparse) (ncol, nrow, anz, TRUE, TRUE, stype,
1096 :     A->xtype, Common) ;
1097 :     if (Common->status < CHOLMOD_OK)
1098 :     {
1099 :     return (FALSE) ; /* out of memory */
1100 :     }
1101 :    
1102 :     if (stype != 0)
1103 :     {
1104 :     /* F = A', upper or lower triangular part only */
1105 :     CHOLMOD(transpose_sym) (A, 1, NULL, F, Common) ;
1106 :     A->packed = TRUE ;
1107 :     /* A = F' */
1108 :     CHOLMOD(transpose_sym) (F, 1, NULL, A, Common) ;
1109 :     }
1110 :     else
1111 :     {
1112 :     /* F = A' */
1113 :     CHOLMOD(transpose_unsym) (A, 1, NULL, NULL, 0, F, Common) ;
1114 :     A->packed = TRUE ;
1115 :     /* A = F' */
1116 :     CHOLMOD(transpose_unsym) (F, 1, NULL, NULL, 0, A, Common) ;
1117 :     }
1118 :    
1119 :     ASSERT (A->sorted && A->packed) ;
1120 :     ASSERT (CHOLMOD(dump_sparse) (A, "Asorted", Common) >= 0) ;
1121 :    
1122 :     /* ---------------------------------------------------------------------- */
1123 :     /* reduce A in size, if needed. This must succeed. */
1124 :     /* ---------------------------------------------------------------------- */
1125 :    
1126 :     Ap = A->p ;
1127 :     anz = Ap [ncol] ;
1128 :     ASSERT ((size_t) anz <= A->nzmax) ;
1129 :     CHOLMOD(reallocate_sparse) (anz, A, Common) ;
1130 :     ASSERT (Common->status >= CHOLMOD_OK) ;
1131 :    
1132 :     /* ---------------------------------------------------------------------- */
1133 :     /* free workspace */
1134 :     /* ---------------------------------------------------------------------- */
1135 :    
1136 :     CHOLMOD(free_sparse) (&F, Common) ;
1137 :     return (TRUE) ;
1138 :     }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge