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 912 - (view) (download) (as text)
Original Path: pkg/src/CHOLMOD/Core/cholmod_transpose.c

1 : bates 912 /* ========================================================================== */
2 :     /* === Core/cholmod_transpose =============================================== */
3 :     /* ========================================================================== */
4 :    
5 :     /* -----------------------------------------------------------------------------
6 :     * CHOLMOD/Core Module. Version 0.6. Copyright (C) 2005, Univ. of Florida.
7 :     * Author: Timothy A. Davis
8 :     * 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 :     * fset = NULL means ":" in MATLAB. fset is ignored.
46 :     * fset != NULL means f = fset [0..fset-1].
47 :     * 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 :     #include "cholmod_core.h"
117 :     #include "cholmod_internal.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 :     Int nrow, ncol, ok = FALSE, permute, use_fset, Apacked, Fpacked, p, pend,
172 :     i, j, k, Fsorted, nf, jj, jlast ;
173 :    
174 :     /* ---------------------------------------------------------------------- */
175 :     /* check inputs */
176 :     /* ---------------------------------------------------------------------- */
177 :    
178 :     RETURN_IF_NULL_COMMON (FALSE) ;
179 :     RETURN_IF_NULL (A, FALSE) ;
180 :     RETURN_IF_NULL (F, FALSE) ;
181 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
182 :     RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
183 :     if (A->nrow != F->ncol || A->ncol != F->nrow)
184 :     {
185 :     ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
186 :     return (FALSE) ;
187 :     }
188 :     Common->status = CHOLMOD_OK ;
189 :    
190 :     /* ---------------------------------------------------------------------- */
191 :     /* get inputs */
192 :     /* ---------------------------------------------------------------------- */
193 :    
194 :     nf = fsize ;
195 :     use_fset = (fset != NULL) ;
196 :     nrow = A->nrow ;
197 :     ncol = A->ncol ;
198 :    
199 :     Ap = A->p ; /* size A->ncol+1, column pointers of A */
200 :     Ai = A->i ; /* size nz = Ap [A->ncol], row indices of A */
201 :     Anz = A->nz ;
202 :     Apacked = A->packed ;
203 :     ASSERT (IMPLIES (!Apacked, Anz != NULL)) ;
204 :    
205 :     permute = (Perm != NULL) ;
206 :    
207 :     Fp = F->p ; /* size A->nrow+1, row pointers of F */
208 :     Fnz = F->nz ;
209 :     Fpacked = F->packed ;
210 :     ASSERT (IMPLIES (!Fpacked, Fnz != NULL)) ;
211 :    
212 :     nf = (use_fset) ? nf : ncol ;
213 :    
214 :     /* ---------------------------------------------------------------------- */
215 :     /* allocate workspace */
216 :     /* ---------------------------------------------------------------------- */
217 :    
218 :     CHOLMOD(allocate_work) (0, nrow + (fset != NULL) ? ncol : 0, 0, Common) ;
219 :     if (Common->status < CHOLMOD_OK)
220 :     {
221 :     return (FALSE) ; /* out of memory */
222 :     }
223 :    
224 :     Wi = Common->Iwork ; /* size nrow (i/l/l) */
225 :    
226 :     /* ---------------------------------------------------------------------- */
227 :     /* check Perm and fset */
228 :     /* ---------------------------------------------------------------------- */
229 :    
230 :     if (permute)
231 :     {
232 :     for (i = 0 ; i < nrow ; i++)
233 :     {
234 :     Wi [i] = 1 ;
235 :     }
236 :     for (k = 0 ; k < nrow ; k++)
237 :     {
238 :     i = Perm [k] ;
239 :     if (i < 0 || i > nrow || Wi [i] == 0)
240 :     {
241 :     ERROR (CHOLMOD_INVALID, "invalid permutation") ;
242 :     return (FALSE) ;
243 :     }
244 :     Wi [i] = 0 ;
245 :     }
246 :     }
247 :    
248 :     if (use_fset)
249 :     {
250 :     for (j = 0 ; j < ncol ; j++)
251 :     {
252 :     Wi [j] = 1 ;
253 :     }
254 :     for (k = 0 ; k < nf ; k++)
255 :     {
256 :     j = fset [k] ;
257 :     if (j < 0 || j > ncol || Wi [j] == 0)
258 :     {
259 :     ERROR (CHOLMOD_INVALID, "invalid fset") ;
260 :     return (FALSE) ;
261 :     }
262 :     Wi [j] = 0 ;
263 :     }
264 :     }
265 :    
266 :     /* Perm and fset are now valid */
267 :     ASSERT (CHOLMOD(dump_perm) (Perm, nrow, nrow, "Perm", Common)) ;
268 :     ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;
269 :    
270 :     /* ---------------------------------------------------------------------- */
271 :     /* count the entries in each row of A or A(:,f) */
272 :     /* ---------------------------------------------------------------------- */
273 :    
274 :     for (i = 0 ; i < nrow ; i++)
275 :     {
276 :     Wi [i] = 0 ;
277 :     }
278 :    
279 :     jlast = EMPTY ;
280 :     Fsorted = TRUE ;
281 :    
282 :     if (use_fset)
283 :     {
284 :     /* count entries in each row of A(:,f) */
285 :     for (jj = 0 ; jj < nf ; jj++)
286 :     {
287 :     j = fset [jj] ;
288 :     if (j <= jlast)
289 :     {
290 :     Fsorted = FALSE ;
291 :     }
292 :     p = Ap [j] ;
293 :     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
294 :     for ( ; p < pend ; p++)
295 :     {
296 :     Wi [Ai [p]]++ ;
297 :     }
298 :     jlast = j ;
299 :     }
300 :    
301 :     /* save the nz counts if F is unpacked, and recount all of A */
302 :     if (!Fpacked)
303 :     {
304 :     if (permute)
305 :     {
306 :     for (i = 0 ; i < nrow ; i++)
307 :     {
308 :     Fnz [i] = Wi [Perm [i]] ;
309 :     }
310 :     }
311 :     else
312 :     {
313 :     for (i = 0 ; i < nrow ; i++)
314 :     {
315 :     Fnz [i] = Wi [i] ;
316 :     }
317 :     }
318 :     for (i = 0 ; i < nrow ; i++)
319 :     {
320 :     Wi [i] = 0 ;
321 :     }
322 :    
323 :     /* count entries in each row of A */
324 :     for (j = 0 ; j < ncol ; j++)
325 :     {
326 :     p = Ap [j] ;
327 :     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
328 :     for ( ; p < pend ; p++)
329 :     {
330 :     Wi [Ai [p]]++ ;
331 :     }
332 :     }
333 :     }
334 :    
335 :     }
336 :     else
337 :     {
338 :    
339 :     /* count entries in each row of A */
340 :     for (j = 0 ; j < ncol ; j++)
341 :     {
342 :     p = Ap [j] ;
343 :     pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
344 :     for ( ; p < pend ; p++)
345 :     {
346 :     Wi [Ai [p]]++ ;
347 :     }
348 :     }
349 :    
350 :     /* save the nz counts if F is unpacked */
351 :     if (!Fpacked)
352 :     {
353 :     if (permute)
354 :     {
355 :     for (i = 0 ; i < nrow ; i++)
356 :     {
357 :     Fnz [i] = Wi [Perm [i]] ;
358 :     }
359 :     }
360 :     else
361 :     {
362 :     for (i = 0 ; i < nrow ; i++)
363 :     {
364 :     Fnz [i] = Wi [i] ;
365 :     }
366 :     }
367 :     }
368 :     }
369 :    
370 :     /* ---------------------------------------------------------------------- */
371 :     /* compute the row pointers */
372 :     /* ---------------------------------------------------------------------- */
373 :    
374 :     p = 0 ;
375 :     if (permute)
376 :     {
377 :     for (i = 0 ; i < nrow ; i++)
378 :     {
379 :     Fp [i] = p ;
380 :     p += Wi [Perm [i]] ;
381 :     }
382 :     for (i = 0 ; i < nrow ; i++)
383 :     {
384 :     Wi [Perm [i]] = Fp [i] ;
385 :     }
386 :     }
387 :     else
388 :     {
389 :     for (i = 0 ; i < nrow ; i++)
390 :     {
391 :     Fp [i] = p ;
392 :     p += Wi [i] ;
393 :     }
394 :     for (i = 0 ; i < nrow ; i++)
395 :     {
396 :     Wi [i] = Fp [i] ;
397 :     }
398 :     }
399 :     Fp [nrow] = p ;
400 :    
401 :     if (p > (Int) (F->nzmax))
402 :     {
403 :     ERROR (CHOLMOD_INVALID, "F is too small") ;
404 :     return (FALSE) ;
405 :     }
406 :    
407 :     /* ---------------------------------------------------------------------- */
408 :     /* transpose matrix, using template routine */
409 :     /* ---------------------------------------------------------------------- */
410 :    
411 :     if (values == 0 || F->xtype == CHOLMOD_PATTERN)
412 :     {
413 :     ok = p_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
414 :     }
415 :     else if (F->xtype == CHOLMOD_REAL)
416 :     {
417 :     ok = r_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
418 :     }
419 :     else if (F->xtype == CHOLMOD_COMPLEX)
420 :     {
421 :     if (values == 1)
422 :     {
423 :     /* array transpose */
424 :     ok = ct_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
425 :     }
426 :     else
427 :     {
428 :     /* complex conjugate transpose */
429 :     ok = c_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
430 :     }
431 :     }
432 :     else if (F->xtype == CHOLMOD_ZOMPLEX)
433 :     {
434 :     if (values == 1)
435 :     {
436 :     /* array transpose */
437 :     ok = zt_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
438 :     }
439 :     else
440 :     {
441 :     /* complex conjugate transpose */
442 :     ok = z_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
443 :     }
444 :     }
445 :    
446 :     /* ---------------------------------------------------------------------- */
447 :     /* finalize result F */
448 :     /* ---------------------------------------------------------------------- */
449 :    
450 :     if (ok)
451 :     {
452 :     F->sorted = Fsorted ;
453 :     }
454 :     ASSERT (CHOLMOD(dump_sparse) (F, "output F unsym", Common) >= 0) ;
455 :     return (ok) ;
456 :     }
457 :    
458 :    
459 :     /* ========================================================================== */
460 :     /* === cholmod_transpose_sym ================================================ */
461 :     /* ========================================================================== */
462 :    
463 :     /* Compute F = A' or A (p,p)', where A is symmetric and F is already allocated.
464 :     * See cholmod_transpose for a simpler routine.
465 :     *
466 :     * workspace: Iwork (nrow) if Perm NULL, Iwork (2*nrow) if Perm non-NULL.
467 :     */
468 :    
469 :     int CHOLMOD(transpose_sym)
470 :     (
471 :     /* ---- input ---- */
472 :     cholmod_sparse *A, /* matrix to transpose */
473 :     int values, /* 2: complex conj. transpose, 1: array transpose,
474 :     0: do not transpose the numerical values */
475 :     Int *Perm, /* size nrow, if present (can be NULL) */
476 :     /* ---- output --- */
477 :     cholmod_sparse *F, /* F = A' or A(p,p)' */
478 :     /* --------------- */
479 :     cholmod_common *Common
480 :     )
481 :     {
482 :     Int *Ap, *Anz, *Ai, *Fp, *Wi, *Pinv, *Iwork ;
483 :     Int p, pend, packed, upper, permute, jold, n, i, j, k, iold ;
484 :     Int ok = FALSE ;
485 :    
486 :     /* ---------------------------------------------------------------------- */
487 :     /* check inputs */
488 :     /* ---------------------------------------------------------------------- */
489 :    
490 :     RETURN_IF_NULL_COMMON (FALSE) ;
491 :     RETURN_IF_NULL (A, FALSE) ;
492 :     RETURN_IF_NULL (F, FALSE) ;
493 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
494 :     RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
495 :     if (A->nrow != A->ncol || A->stype == 0)
496 :     {
497 :     /* this routine handles square symmetric matrices only */
498 :     ERROR (CHOLMOD_INVALID, "matrix must be symmetric") ;
499 :     return (FALSE) ;
500 :     }
501 :     if (A->nrow != F->ncol || A->ncol != F->nrow)
502 :     {
503 :     ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
504 :     return (FALSE) ;
505 :     }
506 :     Common->status = CHOLMOD_OK ;
507 :    
508 :     /* ---------------------------------------------------------------------- */
509 :     /* get inputs */
510 :     /* ---------------------------------------------------------------------- */
511 :    
512 :     permute = (Perm != NULL) ;
513 :     n = A->nrow ;
514 :     Ap = A->p ; /* size A->ncol+1, column pointers of A */
515 :     Ai = A->i ; /* size nz = Ap [A->ncol], row indices of A */
516 :     Anz = A->nz ;
517 :     packed = A->packed ;
518 :     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
519 :     upper = (A->stype > 0) ;
520 :    
521 :     Fp = F->p ; /* size A->nrow+1, row pointers of F */
522 :    
523 :     /* ---------------------------------------------------------------------- */
524 :     /* allocate workspace */
525 :     /* ---------------------------------------------------------------------- */
526 :    
527 :     CHOLMOD(allocate_work) (0, (Perm != NULL) ? 2*n : n, 0, Common) ;
528 :     if (Common->status < CHOLMOD_OK)
529 :     {
530 :     return (FALSE) ; /* out of memory */
531 :     }
532 :    
533 :     /* ---------------------------------------------------------------------- */
534 :     /* get workspace */
535 :     /* ---------------------------------------------------------------------- */
536 :    
537 :     Iwork = Common->Iwork ;
538 :     Wi = Iwork ; /* size n (i/l/l) */
539 :     Pinv = Iwork + n ; /* size n (i/i/l) , unused if Perm NULL */
540 :    
541 :     /* ---------------------------------------------------------------------- */
542 :     /* check Perm and construct inverse permutation */
543 :     /* ---------------------------------------------------------------------- */
544 :    
545 :     if (permute)
546 :     {
547 :     for (i = 0 ; i < n ; i++)
548 :     {
549 :     Pinv [i] = EMPTY ;
550 :     }
551 :     for (k = 0 ; k < n ; k++)
552 :     {
553 :     i = Perm [k] ;
554 :     if (i < 0 || i > n || Pinv [i] != EMPTY)
555 :     {
556 :     ERROR (CHOLMOD_INVALID, "invalid permutation") ;
557 :     return (FALSE) ;
558 :     }
559 :     Pinv [i] = k ;
560 :     }
561 :     }
562 :    
563 :     /* Perm is now valid */
564 :     ASSERT (CHOLMOD(dump_perm) (Perm, n, n, "Perm", Common)) ;
565 :    
566 :     /* ---------------------------------------------------------------------- */
567 :     /* count the entries in each row of F */
568 :     /* ---------------------------------------------------------------------- */
569 :    
570 :     for (i = 0 ; i < n ; i++)
571 :     {
572 :     Wi [i] = 0 ;
573 :     }
574 :    
575 :     if (packed)
576 :     {
577 :     if (permute)
578 :     {
579 :     if (upper)
580 :     {
581 :     /* packed, permuted, upper */
582 :     for (j = 0 ; j < n ; j++)
583 :     {
584 :     jold = Perm [j] ;
585 :     pend = Ap [jold+1] ;
586 :     for (p = Ap [jold] ; p < pend ; p++)
587 :     {
588 :     iold = Ai [p] ;
589 :     if (iold <= jold)
590 :     {
591 :     i = Pinv [iold] ;
592 :     Wi [MIN (i, j)]++ ;
593 :     }
594 :     }
595 :     }
596 :     }
597 :     else
598 :     {
599 :     /* packed, permuted, lower */
600 :     for (j = 0 ; j < n ; j++)
601 :     {
602 :     jold = Perm [j] ;
603 :     pend = Ap [jold+1] ;
604 :     for (p = Ap [jold] ; p < pend ; p++)
605 :     {
606 :     iold = Ai [p] ;
607 :     if (iold >= jold)
608 :     {
609 :     i = Pinv [iold] ;
610 :     Wi [MAX (i, j)]++ ;
611 :     }
612 :     }
613 :     }
614 :     }
615 :     }
616 :     else
617 :     {
618 :     if (upper)
619 :     {
620 :     /* packed, unpermuted, upper */
621 :     for (j = 0 ; j < n ; j++)
622 :     {
623 :     pend = Ap [j+1] ;
624 :     for (p = Ap [j] ; p < pend ; p++)
625 :     {
626 :     i = Ai [p] ;
627 :     if (i <= j)
628 :     {
629 :     Wi [i]++ ;
630 :     }
631 :     }
632 :     }
633 :     }
634 :     else
635 :     {
636 :     /* packed, unpermuted, lower */
637 :     for (j = 0 ; j < n ; j++)
638 :     {
639 :     pend = Ap [j+1] ;
640 :     for (p = Ap [j] ; p < pend ; p++)
641 :     {
642 :     i = Ai [p] ;
643 :     if (i >= j)
644 :     {
645 :     Wi [i]++ ;
646 :     }
647 :     }
648 :     }
649 :     }
650 :     }
651 :     }
652 :     else
653 :     {
654 :     if (permute)
655 :     {
656 :     if (upper)
657 :     {
658 :     /* unpacked, permuted, upper */
659 :     for (j = 0 ; j < n ; j++)
660 :     {
661 :     jold = Perm [j] ;
662 :     p = Ap [jold] ;
663 :     pend = p + Anz [jold] ;
664 :     for ( ; p < pend ; p++)
665 :     {
666 :     iold = Ai [p] ;
667 :     if (iold <= jold)
668 :     {
669 :     i = Pinv [iold] ;
670 :     Wi [MIN (i, j)]++ ;
671 :     }
672 :     }
673 :     }
674 :     }
675 :     else
676 :     {
677 :     /* unpacked, permuted, lower */
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 [MAX (i, j)]++ ;
690 :     }
691 :     }
692 :     }
693 :     }
694 :     }
695 :     else
696 :     {
697 :     if (upper)
698 :     {
699 :     /* unpacked, unpermuted, upper */
700 :     for (j = 0 ; j < n ; j++)
701 :     {
702 :     p = Ap [j] ;
703 :     pend = p + Anz [j] ;
704 :     for ( ; p < pend ; p++)
705 :     {
706 :     i = Ai [p] ;
707 :     if (i <= j)
708 :     {
709 :     Wi [i]++ ;
710 :     }
711 :     }
712 :     }
713 :     }
714 :     else
715 :     {
716 :     /* unpacked, unpermuted, lower */
717 :     for (j = 0 ; j < n ; j++)
718 :     {
719 :     p = Ap [j] ;
720 :     pend = p + Anz [j] ;
721 :     for ( ; p < pend ; p++)
722 :     {
723 :     i = Ai [p] ;
724 :     if (i >= j)
725 :     {
726 :     Wi [i]++ ;
727 :     }
728 :     }
729 :     }
730 :     }
731 :     }
732 :     }
733 :    
734 :     /* ---------------------------------------------------------------------- */
735 :     /* compute the row pointers */
736 :     /* ---------------------------------------------------------------------- */
737 :    
738 :     p = 0 ;
739 :     for (i = 0 ; i < n ; i++)
740 :     {
741 :     Fp [i] = p ;
742 :     p += Wi [i] ;
743 :     }
744 :     Fp [n] = p ;
745 :     for (i = 0 ; i < n ; i++)
746 :     {
747 :     Wi [i] = Fp [i] ;
748 :     }
749 :    
750 :     if (p > (Int) (F->nzmax))
751 :     {
752 :     ERROR (CHOLMOD_INVALID, "F is too small") ;
753 :     return (FALSE) ;
754 :     }
755 :    
756 :     /* ---------------------------------------------------------------------- */
757 :     /* transpose matrix, using template routine */
758 :     /* ---------------------------------------------------------------------- */
759 :    
760 :     if (values == 0 || F->xtype == CHOLMOD_PATTERN)
761 :     {
762 :     PRINT2 (("\n:::: p_transpose_sym Perm %p\n", Perm)) ;
763 :     ok = p_cholmod_transpose_sym (A, Perm, F, Common) ;
764 :     }
765 :     else if (F->xtype == CHOLMOD_REAL)
766 :     {
767 :     PRINT2 (("\n:::: r_transpose_sym Perm %p\n", Perm)) ;
768 :     ok = r_cholmod_transpose_sym (A, Perm, F, Common) ;
769 :     }
770 :     else if (F->xtype == CHOLMOD_COMPLEX)
771 :     {
772 :     if (values == 1)
773 :     {
774 :     /* array transpose */
775 :     PRINT2 (("\n:::: ct_transpose_sym Perm %p\n", Perm)) ;
776 :     ok = ct_cholmod_transpose_sym (A, Perm, F, Common) ;
777 :     }
778 :     else
779 :     {
780 :     /* complex conjugate transpose */
781 :     PRINT2 (("\n:::: c_transpose_sym Perm %p\n", Perm)) ;
782 :     ok = c_cholmod_transpose_sym (A, Perm, F, Common) ;
783 :     }
784 :     }
785 :     else if (F->xtype == CHOLMOD_ZOMPLEX)
786 :     {
787 :     if (values == 1)
788 :     {
789 :     /* array transpose */
790 :     PRINT2 (("\n:::: zt_transpose_sym Perm %p\n", Perm)) ;
791 :     ok = zt_cholmod_transpose_sym (A, Perm, F, Common) ;
792 :     }
793 :     else
794 :     {
795 :     /* complex conjugate transpose */
796 :     PRINT2 (("\n:::: z_transpose_sym Perm %p\n", Perm)) ;
797 :     ok = z_cholmod_transpose_sym (A, Perm, F, Common) ;
798 :     }
799 :     }
800 :    
801 :     /* ---------------------------------------------------------------------- */
802 :     /* finalize result F */
803 :     /* ---------------------------------------------------------------------- */
804 :    
805 :     /* F is sorted if there is no permutation vector */
806 :     if (ok)
807 :     {
808 :     F->sorted = !permute ;
809 :     F->packed = TRUE ;
810 :     F->stype = - SIGN (A->stype) ; /* flip the stype */
811 :     ASSERT (CHOLMOD(dump_sparse) (F, "output F sym", Common) >= 0) ;
812 :     }
813 :     return (ok) ;
814 :     }
815 :    
816 :    
817 :     /* ========================================================================== */
818 :     /* === cholmod_transpose ==================================================== */
819 :     /* ========================================================================== */
820 :    
821 :     /* Returns A'. See also cholmod_ptranspose below. */
822 :    
823 :     cholmod_sparse *CHOLMOD(transpose)
824 :     (
825 :     /* ---- input ---- */
826 :     cholmod_sparse *A, /* matrix to transpose */
827 :     int values, /* 2: complex conj. transpose, 1: array transpose,
828 :     0: do not transpose the numerical values
829 :     (returns its result as CHOLMOD_PATTERN) */
830 :     /* --------------- */
831 :     cholmod_common *Common
832 :     )
833 :     {
834 :     return (CHOLMOD(ptranspose) (A, values, NULL, NULL, 0, Common)) ;
835 :     }
836 :    
837 :    
838 :     /* ========================================================================== */
839 :     /* === cholmod_ptranspose =================================================== */
840 :     /* ========================================================================== */
841 :    
842 :     /* Return A' or A(p,p)' if A is symmetric. Return A', A(:,f)', or A(p,f)' if
843 :     * A is unsymmetric.
844 :     *
845 :     * workspace:
846 :     * Iwork (MAX (nrow,ncol)) if unsymmetric and fset is non-NULL
847 :     * Iwork (nrow) if unsymmetric and fset is NULL
848 :     * Iwork (2*nrow) if symmetric and Perm is non-NULL.
849 :     * Iwork (nrow) if symmetric and Perm is NULL.
850 :     *
851 :     * A simple worst-case upper bound on the workspace is nrow+ncol.
852 :     */
853 :    
854 :     cholmod_sparse *CHOLMOD(ptranspose)
855 :     (
856 :     /* ---- input ---- */
857 :     cholmod_sparse *A, /* matrix to transpose */
858 :     int values, /* 2: complex conj. transpose, 1: array transpose,
859 :     0: do not transpose the numerical values */
860 :     Int *Perm, /* if non-NULL, F = A(p,f) or A(p,p) */
861 :     Int *fset, /* subset of 0:(A->ncol)-1 */
862 :     size_t fsize, /* size of fset */
863 :     /* --------------- */
864 :     cholmod_common *Common
865 :     )
866 :     {
867 :     Int *Ap, *Anz ;
868 :     cholmod_sparse *F ;
869 :     Int nrow, ncol, use_fset, j, jj, fnz, packed, stype, nf, ineed, ok,
870 :     xtype ;
871 :    
872 :     nf = fsize ;
873 :    
874 :     /* ---------------------------------------------------------------------- */
875 :     /* check inputs */
876 :     /* ---------------------------------------------------------------------- */
877 :    
878 :     RETURN_IF_NULL_COMMON (NULL) ;
879 :     RETURN_IF_NULL (A, FALSE) ;
880 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
881 :     stype = A->stype ;
882 :     Common->status = CHOLMOD_OK ;
883 :    
884 :     /* ---------------------------------------------------------------------- */
885 :     /* allocate workspace */
886 :     /* ---------------------------------------------------------------------- */
887 :    
888 :     nrow = A->nrow ;
889 :     ncol = A->ncol ;
890 :    
891 :     if (stype != 0)
892 :     {
893 :     use_fset = FALSE ;
894 :     if (Perm != NULL)
895 :     {
896 :     ineed = 2*nrow ;
897 :     }
898 :     else
899 :     {
900 :     ineed = nrow ;
901 :     }
902 :     }
903 :     else
904 :     {
905 :     use_fset = (fset != NULL) ;
906 :     if (use_fset)
907 :     {
908 :     ineed = MAX (nrow, ncol) ;
909 :     }
910 :     else
911 :     {
912 :     ineed = nrow ;
913 :     }
914 :     }
915 :    
916 :     CHOLMOD(allocate_work) (0, ineed, 0, Common) ;
917 :     if (Common->status < CHOLMOD_OK)
918 :     {
919 :     return (NULL) ; /* out of memory */
920 :     }
921 :    
922 :     /* ---------------------------------------------------------------------- */
923 :     /* get inputs */
924 :     /* ---------------------------------------------------------------------- */
925 :    
926 :     Ap = A->p ;
927 :     Anz = A->nz ;
928 :     packed = A->packed ;
929 :     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
930 :     xtype = values ? A->xtype : CHOLMOD_PATTERN ;
931 :    
932 :     /* ---------------------------------------------------------------------- */
933 :     /* allocate F */
934 :     /* ---------------------------------------------------------------------- */
935 :    
936 :     /* determine # of nonzeros in F */
937 :     if (stype != 0)
938 :     {
939 :     /* F=A' or F=A(p,p)', fset is ignored */
940 :     fnz = CHOLMOD(nnz) (A, Common) ;
941 :     }
942 :     else
943 :     {
944 :     nf = (use_fset) ? nf : ncol ;
945 :     if (use_fset)
946 :     {
947 :     fnz = 0 ;
948 :     /* F=A(:,f)' or F=A(p,f)' */
949 :     for (jj = 0 ; jj < nf ; jj++)
950 :     {
951 :     /* The fset is not yet checked; it will be thoroughly checked
952 :     * in cholmod_transpose_unsym. For now, just make sure we don't
953 :     * access Ap and Anz out of bounds. */
954 :     j = fset [jj] ;
955 :     if (j >= 0 && j < ncol)
956 :     {
957 :     fnz += packed ? (Ap [j+1] - Ap [j]) : MAX (0, Anz [j]) ;
958 :     }
959 :     }
960 :     }
961 :     else
962 :     {
963 :     /* F=A' or F=A(p,:)' */
964 :     fnz = CHOLMOD(nnz) (A, Common) ;
965 :     }
966 :     }
967 :    
968 :     /* F is ncol-by-nrow, fnz nonzeros, sorted unless f is present and unsorted,
969 :     * packed, of opposite stype as A, and with/without numerical values */
970 :     F = CHOLMOD(allocate_sparse) (ncol, nrow, fnz, TRUE, TRUE, -SIGN(stype),
971 :     xtype, Common) ;
972 :     if (Common->status < CHOLMOD_OK)
973 :     {
974 :     return (NULL) ; /* out of memory */
975 :     }
976 :    
977 :     /* ---------------------------------------------------------------------- */
978 :     /* transpose and optionally permute the matrix A */
979 :     /* ---------------------------------------------------------------------- */
980 :    
981 :     if (stype != 0)
982 :     {
983 :     /* F = A (p,p)', using upper or lower triangular part of A only */
984 :     ok = CHOLMOD(transpose_sym) (A, values, Perm, F, Common) ;
985 :     }
986 :     else
987 :     {
988 :     /* F = A (p,f)' */
989 :     ok = CHOLMOD(transpose_unsym) (A, values, Perm, fset, nf, F, Common) ;
990 :     }
991 :    
992 :     /* ---------------------------------------------------------------------- */
993 :     /* return the matrix F, or NULL if an error occured */
994 :     /* ---------------------------------------------------------------------- */
995 :    
996 :     if (!ok)
997 :     {
998 :     CHOLMOD(free_sparse) (&F, Common) ;
999 :     }
1000 :     return (F) ;
1001 :     }
1002 :    
1003 :    
1004 :     /* ========================================================================== */
1005 :     /* === cholmod_sort ========================================================= */
1006 :     /* ========================================================================== */
1007 :    
1008 :     /* Sort the columns of A, in place. Returns A in packed form, even if it
1009 :     * starts as unpacked. Removes entries in the ignored part of a symmetric
1010 :     * matrix.
1011 :     *
1012 :     * workspace: Iwork (max (nrow,ncol)). Allocates additional workspace for a
1013 :     * temporary copy of A'.
1014 :     */
1015 :    
1016 :     int CHOLMOD(sort)
1017 :     (
1018 :     /* ---- in/out --- */
1019 :     cholmod_sparse *A, /* matrix to sort */
1020 :     /* --------------- */
1021 :     cholmod_common *Common
1022 :     )
1023 :     {
1024 :     Int *Ap ;
1025 :     cholmod_sparse *F ;
1026 :     Int anz, ncol, nrow, stype ;
1027 :    
1028 :     /* ---------------------------------------------------------------------- */
1029 :     /* check inputs */
1030 :     /* ---------------------------------------------------------------------- */
1031 :    
1032 :     RETURN_IF_NULL_COMMON (FALSE) ;
1033 :     RETURN_IF_NULL (A, FALSE) ;
1034 :     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
1035 :     Common->status = CHOLMOD_OK ;
1036 :     nrow = A->nrow ;
1037 :     if (nrow <= 1)
1038 :     {
1039 :     /* a 1-by-n sparse matrix must be sorted */
1040 :     A->sorted = TRUE ;
1041 :     return (TRUE) ;
1042 :     }
1043 :    
1044 :     /* ---------------------------------------------------------------------- */
1045 :     /* allocate workspace */
1046 :     /* ---------------------------------------------------------------------- */
1047 :    
1048 :     ncol = A->ncol ;
1049 :     CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
1050 :     if (Common->status < CHOLMOD_OK)
1051 :     {
1052 :     return (FALSE) ; /* out of memory */
1053 :     }
1054 :    
1055 :     /* ---------------------------------------------------------------------- */
1056 :     /* get inputs */
1057 :     /* ---------------------------------------------------------------------- */
1058 :    
1059 :     anz = CHOLMOD(nnz) (A, Common) ;
1060 :     stype = A->stype ;
1061 :    
1062 :     /* ---------------------------------------------------------------------- */
1063 :     /* sort the columns of the matrix */
1064 :     /* ---------------------------------------------------------------------- */
1065 :    
1066 :     /* allocate workspace for transpose: ncol-by-nrow, same # of nonzeros as A,
1067 :     * sorted, packed, same stype as A, and of the same numeric type as A. */
1068 :     F = CHOLMOD(allocate_sparse) (ncol, nrow, anz, TRUE, TRUE, stype,
1069 :     A->xtype, Common) ;
1070 :     if (Common->status < CHOLMOD_OK)
1071 :     {
1072 :     return (FALSE) ; /* out of memory */
1073 :     }
1074 :    
1075 :     if (stype != 0)
1076 :     {
1077 :     /* F = A', upper or lower triangular part only */
1078 :     CHOLMOD(transpose_sym) (A, 1, NULL, F, Common) ;
1079 :     A->packed = TRUE ;
1080 :     /* A = F' */
1081 :     CHOLMOD(transpose_sym) (F, 1, NULL, A, Common) ;
1082 :     }
1083 :     else
1084 :     {
1085 :     /* F = A' */
1086 :     CHOLMOD(transpose_unsym) (A, 1, NULL, NULL, 0, F, Common) ;
1087 :     A->packed = TRUE ;
1088 :     /* A = F' */
1089 :     CHOLMOD(transpose_unsym) (F, 1, NULL, NULL, 0, A, Common) ;
1090 :     }
1091 :    
1092 :     ASSERT (A->sorted && A->packed) ;
1093 :     ASSERT (CHOLMOD(dump_sparse) (A, "Asorted", Common) >= 0) ;
1094 :    
1095 :     /* ---------------------------------------------------------------------- */
1096 :     /* reduce A in size, if needed. This must succeed. */
1097 :     /* ---------------------------------------------------------------------- */
1098 :    
1099 :     Ap = A->p ;
1100 :     anz = Ap [ncol] ;
1101 :     ASSERT ((size_t) anz <= A->nzmax) ;
1102 :     CHOLMOD(reallocate_sparse) (anz, A, Common) ;
1103 :     ASSERT (Common->status >= CHOLMOD_OK) ;
1104 :    
1105 :     /* ---------------------------------------------------------------------- */
1106 :     /* free workspace */
1107 :     /* ---------------------------------------------------------------------- */
1108 :    
1109 :     CHOLMOD(free_sparse) (&F, Common) ;
1110 :     return (TRUE) ;
1111 :     }

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