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

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