1 /* ========================================================================== */
2 /* === Core/cholmod_triplet ================================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Core Module. Copyright (C) 2005-2006,
7 * Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* Core utility routines for the cholmod_triplet object:
11 *
12 * A sparse matrix held in triplet form is the simplest one for a user to
13 * create. It consists of a list of nz entries in arbitrary order, held in
14 * three arrays: i, j, and x, each of length nk. The kth entry is in row i[k],
15 * column j[k], with value x[k]. There may be duplicate values; if A(i,j)
16 * appears more than once, its value is the sum of the entries with those row
17 * and column indices.
18 *
19 * Primary routines:
20 * -----------------
21 * cholmod_allocate_triplet allocate a triplet matrix
22 * cholmod_free_triplet free a triplet matrix
23 *
24 * Secondary routines:
25 * -------------------
26 * cholmod_reallocate_triplet reallocate a triplet matrix
27 * cholmod_sparse_to_triplet create a triplet matrix copy of a sparse matrix
28 * cholmod_triplet_to_sparse create a sparse matrix copy of a triplet matrix
29 * cholmod_copy_triplet create a copy of a triplet matrix
30 *
31 * The relationship between an m-by-n cholmod_sparse matrix A and a
32 * cholmod_triplet matrix (i, j, and x) is identical to how they are used in
33 * the MATLAB "sparse" and "find" functions:
34 *
35 * [i j x] = find (A)
36 * [m n] = size (A)
37 * A = sparse (i,j,x,m,n)
38 *
39 * with the exception that the cholmod_sparse matrix may be "unpacked", may
40 * have either sorted or unsorted columns (depending on the option selected),
41 * and may be symmetric with just the upper or lower triangular part stored.
42 * Likewise, the cholmod_triplet matrix may contain just the entries in the
43 * upper or lower triangular part of a symmetric matrix.
44 *
45 * MATLAB sparse matrices are always "packed", always have sorted columns,
46 * and always store both parts of a symmetric matrix. In some cases, MATLAB
47 * behaves like CHOLMOD by ignoring entries in the upper or lower triangular
48 * part of a matrix that is otherwise assumed to be symmetric (such as the
49 * input to chol). In CHOLMOD, that option is a characteristic of the object.
50 * In MATLAB, that option is based on how a matrix is used as the input to
51 * a function.
52 *
53 * The triplet matrix is provided to give the user a simple way of constructing
54 * a sparse matrix. There are very few operations supported for triplet
55 * matrices. The assumption is that they will be converted to cholmod_sparse
56 * matrix form first.
57 *
58 * Adding two triplet matrices simply involves concatenating the contents of
59 * the three arrays (i, j, and x). To permute a triplet matrix, just replace
60 * the row and column indices with their permuted values. For example, if
61 * P is a permutation vector, then P [k] = j means row/column j is the kth
62 * row/column in C=P*A*P'. In MATLAB notation, C=A(p,p). If Pinv is an array
63 * of size n and T is the triplet form of A, then:
64 *
65 * Ti = T->i ;
66 * Tj = T->j ;
67 * for (k = 0 ; k < n ; k++) Pinv [P [k]] = k ;
68 * for (k = 0 ; k < nz ; k++) Ti [k] = Pinv [Ti [k]] ;
69 * for (k = 0 ; k < nz ; k++) Tj [k] = Pinv [Tj [k]] ;
70 *
71 * overwrites T with the triplet form of C=P*A*P'. The conversion
72 *
73 * C = cholmod_triplet_to_sparse (T, 0, &Common) ;
74 *
75 * will then return the matrix C = P*A*P'.
76 *
77 * Note that T->stype > 0 means that entries in the lower triangular part of
78 * T are transposed into the upper triangular part when T is converted to
79 * sparse matrix (cholmod_sparse) form with cholmod_triplet_to_sparse. The
80 * opposite is true for T->stype < 0.
81 *
82 * Since the triplet matrix T is so simple to generate, it's quite easy
83 * to remove entries that you do not want, prior to converting T to the
84 * cholmod_sparse form. So if you include these entries in T, CHOLMOD
85 * assumes that there must be a reason (such as the one above). Thus,
86 * no entry in a triplet matrix is ever ignored.
87 *
88 * Other operations, such as extacting a submatrix, horizontal and vertical
89 * concatenation, multiply a triplet matrix times a dense matrix, are also
90 * simple. Multiplying two triplet matrices is not trivial; the simplest
91 * method is to convert them to cholmod_sparse matrices first.
92 *
93 * Supports all xtypes (pattern, real, complex, and zomplex).
94 */
95
96 #include "cholmod_internal.h"
97 #include "cholmod_core.h"
98
99
100 /* ========================================================================== */
101 /* === TEMPLATE ============================================================= */
102 /* ========================================================================== */
103
104 #define PATTERN
105 #include "t_cholmod_triplet.c"
106 #define REAL
107 #include "t_cholmod_triplet.c"
108 #define COMPLEX
109 #include "t_cholmod_triplet.c"
110 #define ZOMPLEX
111 #include "t_cholmod_triplet.c"
112
113
114 /* ========================================================================== */
115 /* === cholmod_allocate_triplet ============================================= */
116 /* ========================================================================== */
117
118 /* allocate space for a triplet matrix
119 *
120 * workspace: none
121 */
122
CHOLMOD(allocate_triplet)123 cholmod_triplet *CHOLMOD(allocate_triplet)
124 (
125 /* ---- input ---- */
126 size_t nrow, /* # of rows of T */
127 size_t ncol, /* # of columns of T */
128 size_t nzmax, /* max # of nonzeros of T */
129 int stype, /* stype of T */
130 int xtype, /* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
131 /* --------------- */
132 cholmod_common *Common
133 )
134 {
135 cholmod_triplet *T ;
136 size_t nzmax0 ;
137 int ok = TRUE ;
138
139 /* ---------------------------------------------------------------------- */
140 /* get inputs */
141 /* ---------------------------------------------------------------------- */
142
143 RETURN_IF_NULL_COMMON (NULL) ;
144 if (xtype < CHOLMOD_PATTERN || xtype > CHOLMOD_ZOMPLEX)
145 {
146 ERROR (CHOLMOD_INVALID, "xtype invalid") ;
147 return (NULL) ;
148 }
149 /* ensure the dimensions do not cause integer overflow */
150 (void) CHOLMOD(add_size_t) (ncol, 2, &ok) ;
151 if (!ok || nrow > Int_max || ncol > Int_max || nzmax > Int_max)
152 {
153 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
154 return (NULL) ;
155 }
156
157 Common->status = CHOLMOD_OK ;
158
159 /* ---------------------------------------------------------------------- */
160 /* allocate header */
161 /* ---------------------------------------------------------------------- */
162
163 T = CHOLMOD(malloc) (sizeof (cholmod_triplet), 1, Common) ;
164 if (Common->status < CHOLMOD_OK)
165 {
166 return (NULL) ; /* out of memory */
167 }
168
169 PRINT1 (("cholmod_allocate_triplet %d-by-%d nzmax %d xtype %d\n",
170 nrow, ncol, nzmax, xtype)) ;
171
172 nzmax = MAX (1, nzmax) ;
173
174 T->nrow = nrow ;
175 T->ncol = ncol ;
176 T->nzmax = nzmax ;
177 T->nnz = 0 ;
178 T->stype = stype ;
179 T->itype = ITYPE ;
180 T->xtype = xtype ;
181 T->dtype = DTYPE ;
182
183 T->j = NULL ;
184 T->i = NULL ;
185 T->x = NULL ;
186 T->z = NULL ;
187
188 /* ---------------------------------------------------------------------- */
189 /* allocate the matrix itself */
190 /* ---------------------------------------------------------------------- */
191
192 nzmax0 = 0 ;
193 CHOLMOD(realloc_multiple) (nzmax, 2, xtype, &(T->i), &(T->j),
194 &(T->x), &(T->z), &nzmax0, Common) ;
195
196 if (Common->status < CHOLMOD_OK)
197 {
198 CHOLMOD(free_triplet) (&T, Common) ;
199 return (NULL) ; /* out of memory */
200 }
201
202 return (T) ;
203 }
204
205
206 /* ========================================================================== */
207 /* === cholmod_free_triplet ================================================= */
208 /* ========================================================================== */
209
210 /* free a triplet matrix
211 *
212 * workspace: none
213 */
214
CHOLMOD(free_triplet)215 int CHOLMOD(free_triplet)
216 (
217 /* ---- in/out --- */
218 cholmod_triplet **THandle, /* matrix to deallocate, NULL on output */
219 /* --------------- */
220 cholmod_common *Common
221 )
222 {
223 Int nz ;
224 cholmod_triplet *T ;
225
226 RETURN_IF_NULL_COMMON (FALSE) ;
227
228 if (THandle == NULL)
229 {
230 /* nothing to do */
231 return (TRUE) ;
232 }
233 T = *THandle ;
234 if (T == NULL)
235 {
236 /* nothing to do */
237 return (TRUE) ;
238 }
239 nz = T->nzmax ;
240 T->j = CHOLMOD(free) (nz, sizeof (Int), T->j, Common) ;
241 T->i = CHOLMOD(free) (nz, sizeof (Int), T->i, Common) ;
242 if (T->xtype == CHOLMOD_REAL)
243 {
244 T->x = CHOLMOD(free) (nz, sizeof (double), T->x, Common) ;
245 }
246 else if (T->xtype == CHOLMOD_COMPLEX)
247 {
248 T->x = CHOLMOD(free) (nz, 2*sizeof (double), T->x, Common) ;
249 }
250 else if (T->xtype == CHOLMOD_ZOMPLEX)
251 {
252 T->x = CHOLMOD(free) (nz, sizeof (double), T->x, Common) ;
253 T->z = CHOLMOD(free) (nz, sizeof (double), T->z, Common) ;
254 }
255 *THandle = CHOLMOD(free) (1, sizeof (cholmod_triplet), (*THandle), Common) ;
256 return (TRUE) ;
257 }
258
259
260 /* ========================================================================== */
261 /* === cholmod_reallocate_triplet =========================================== */
262 /* ========================================================================== */
263
264 /* Change the size of T->i, T->j, and T->x, or allocate them if their current
265 * size is zero. T->x is not modified if T->xtype is CHOLMOD_PATTERN.
266 *
267 * workspace: none
268 */
269
CHOLMOD(reallocate_triplet)270 int CHOLMOD(reallocate_triplet)
271 (
272 /* ---- input ---- */
273 size_t nznew, /* new # of entries in T */
274 /* ---- in/out --- */
275 cholmod_triplet *T, /* triplet matrix to modify */
276 /* --------------- */
277 cholmod_common *Common
278 )
279 {
280
281 /* ---------------------------------------------------------------------- */
282 /* get inputs */
283 /* ---------------------------------------------------------------------- */
284
285 RETURN_IF_NULL_COMMON (FALSE) ;
286 RETURN_IF_NULL (T, FALSE) ;
287 RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
288 Common->status = CHOLMOD_OK ;
289 PRINT1 (("realloc triplet %d to %d, xtype: %d\n",
290 T->nzmax, nznew, T->xtype)) ;
291
292 /* ---------------------------------------------------------------------- */
293 /* resize the matrix */
294 /* ---------------------------------------------------------------------- */
295
296 CHOLMOD(realloc_multiple) (MAX (1,nznew), 2, T->xtype, &(T->i), &(T->j),
297 &(T->x), &(T->z), &(T->nzmax), Common) ;
298
299 return (Common->status == CHOLMOD_OK) ;
300 }
301
302
303 /* ========================================================================== */
304 /* === cholmod_triplet_to_sparse ============================================ */
305 /* ========================================================================== */
306
307 /* Convert a set of triplets into a cholmod_sparse matrix. In MATLAB notation,
308 * for unsymmetric matrices:
309 *
310 * A = sparse (Ti, Tj, Tx, nrow, ncol, nzmax) ;
311 *
312 * For the symmetric upper case:
313 *
314 * A = sparse (min(Ti,Tj), max(Ti,Tj), Tx, nrow, ncol, nzmax) ;
315 *
316 * For the symmetric lower case:
317 *
318 * A = sparse (max(Ti,Tj), min(Ti,Tj), Tx, nrow, ncol, nzmax) ;
319 *
320 * If Tx is NULL, then A->x is not allocated, and only the pattern of A is
321 * computed. A is returned in packed form, and can be of any stype
322 * (upper/lower/unsymmetric). It has enough space to hold the values in T,
323 * or nzmax, whichever is larger.
324 *
325 * workspace: Iwork (max (nrow,ncol))
326 * allocates a temporary copy of its output matrix.
327 *
328 * The resulting sparse matrix has the same xtype as the input triplet matrix.
329 */
330
CHOLMOD(triplet_to_sparse)331 cholmod_sparse *CHOLMOD(triplet_to_sparse)
332 (
333 /* ---- input ---- */
334 cholmod_triplet *T, /* matrix to copy */
335 size_t nzmax, /* allocate at least this much space in output matrix */
336 /* --------------- */
337 cholmod_common *Common
338 )
339 {
340 cholmod_sparse *R, *A = NULL ;
341 Int *Wj, *Rp, *Ri, *Rnz, *Ti, *Tj ;
342 Int i, j, p, k, stype, nrow, ncol, nz, ok ;
343 size_t anz = 0 ;
344
345 /* ---------------------------------------------------------------------- */
346 /* check inputs */
347 /* ---------------------------------------------------------------------- */
348
349 RETURN_IF_NULL_COMMON (NULL) ;
350 RETURN_IF_NULL (T, NULL) ;
351 Ti = T->i ;
352 Tj = T->j ;
353 RETURN_IF_NULL (Ti, NULL) ;
354 RETURN_IF_NULL (Tj, NULL) ;
355 RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
356 stype = SIGN (T->stype) ;
357 if (stype && T->nrow != T->ncol)
358 {
359 /* inputs invalid */
360 ERROR (CHOLMOD_INVALID, "matrix invalid") ;
361 return (NULL) ;
362 }
363 Common->status = CHOLMOD_OK ;
364 DEBUG (CHOLMOD(dump_triplet) (T, "T", Common)) ;
365
366 /* ---------------------------------------------------------------------- */
367 /* get inputs */
368 /* ---------------------------------------------------------------------- */
369
370 nrow = T->nrow ;
371 ncol = T->ncol ;
372 nz = T->nnz ;
373
374 /* ---------------------------------------------------------------------- */
375 /* allocate workspace */
376 /* ---------------------------------------------------------------------- */
377
378 CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
379 if (Common->status < CHOLMOD_OK)
380 {
381 return (NULL) ; /* out of memory */
382 }
383
384 /* ---------------------------------------------------------------------- */
385 /* allocate temporary matrix R */
386 /* ---------------------------------------------------------------------- */
387
388 R = CHOLMOD(allocate_sparse) (ncol, nrow, nz, FALSE, FALSE, -stype,
389 T->xtype, Common) ;
390
391 if (Common->status < CHOLMOD_OK)
392 {
393 return (NULL) ; /* out of memory */
394 }
395
396 Rp = R->p ;
397 Ri = R->i ;
398 Rnz = R->nz ;
399
400 /* ---------------------------------------------------------------------- */
401 /* count the entries in each row of A (also counting duplicates) */
402 /* ---------------------------------------------------------------------- */
403
404 for (i = 0 ; i < nrow ; i++)
405 {
406 Rnz [i] = 0 ;
407 }
408
409 if (stype > 0)
410 {
411 for (k = 0 ; k < nz ; k++)
412 {
413 i = Ti [k] ;
414 j = Tj [k] ;
415 if (i < 0 || i >= nrow || j < 0 || j >= ncol)
416 {
417 ERROR (CHOLMOD_INVALID, "index out of range") ;
418 break ;
419 }
420 /* A will be symmetric with just the upper triangular part stored.
421 * Create a matrix R that is lower triangular. Entries in the
422 * upper part of R are transposed to the lower part. */
423 Rnz [MIN (i,j)]++ ;
424 }
425 }
426 else if (stype < 0)
427 {
428 for (k = 0 ; k < nz ; k++)
429 {
430 i = Ti [k] ;
431 j = Tj [k] ;
432 if (i < 0 || i >= nrow || j < 0 || j >= ncol)
433 {
434 ERROR (CHOLMOD_INVALID, "index out of range") ;
435 break ;
436 }
437 /* A will be symmetric with just the lower triangular part stored.
438 * Create a matrix R that is upper triangular. Entries in the
439 * lower part of R are transposed to the upper part. */
440 Rnz [MAX (i,j)]++ ;
441 }
442 }
443 else
444 {
445 for (k = 0 ; k < nz ; k++)
446 {
447 i = Ti [k] ;
448 j = Tj [k] ;
449 if (i < 0 || i >= nrow || j < 0 || j >= ncol)
450 {
451 ERROR (CHOLMOD_INVALID, "index out of range") ;
452 break ;
453 }
454 /* constructing an unsymmetric matrix */
455 Rnz [i]++ ;
456 }
457 }
458
459 if (Common->status < CHOLMOD_OK)
460 {
461 /* triplet matrix is invalid */
462 CHOLMOD(free_sparse) (&R, Common) ;
463 return (NULL) ;
464 }
465
466 /* ---------------------------------------------------------------------- */
467 /* construct the row pointers */
468 /* ---------------------------------------------------------------------- */
469
470 p = 0 ;
471 for (i = 0 ; i < nrow ; i++)
472 {
473 Rp [i] = p ;
474 p += Rnz [i] ;
475 }
476 Rp [nrow] = p ;
477
478 /* use Wj (i/l/l) as temporary row pointers */
479 Wj = Common->Iwork ; /* size MAX (nrow,ncol) FUTURE WORK: (i/l/l) */
480 for (i = 0 ; i < nrow ; i++)
481 {
482 Wj [i] = Rp [i] ;
483 }
484
485 /* ---------------------------------------------------------------------- */
486 /* construct triplet matrix, using template routine */
487 /* ---------------------------------------------------------------------- */
488
489 switch (T->xtype)
490 {
491 case CHOLMOD_PATTERN:
492 anz = p_cholmod_triplet_to_sparse (T, R, Common) ;
493 break ;
494
495 case CHOLMOD_REAL:
496 anz = r_cholmod_triplet_to_sparse (T, R, Common) ;
497 break ;
498
499 case CHOLMOD_COMPLEX:
500 anz = c_cholmod_triplet_to_sparse (T, R, Common) ;
501 break ;
502
503 case CHOLMOD_ZOMPLEX:
504 anz = z_cholmod_triplet_to_sparse (T, R, Common) ;
505 break ;
506 }
507
508 /* ---------------------------------------------------------------------- */
509 /* A = R' (array transpose, not complex conjugate transpose) */
510 /* ---------------------------------------------------------------------- */
511
512 /* workspace: Iwork (R->nrow), which is A->ncol */
513
514 ASSERT (CHOLMOD(dump_sparse) (R, "R", Common) >= 0) ;
515
516 A = CHOLMOD(allocate_sparse) (nrow, ncol, MAX (anz, nzmax), TRUE, TRUE,
517 stype, T->xtype, Common) ;
518
519 if (stype)
520 {
521 ok = CHOLMOD(transpose_sym) (R, 1, NULL, A, Common) ;
522 }
523 else
524 {
525 ok = CHOLMOD(transpose_unsym) (R, 1, NULL, NULL, 0, A, Common) ;
526 }
527
528 CHOLMOD(free_sparse) (&R, Common) ;
529 if (Common->status < CHOLMOD_OK)
530 {
531 CHOLMOD(free_sparse) (&A, Common) ;
532 }
533
534 /* ---------------------------------------------------------------------- */
535 /* return result */
536 /* ---------------------------------------------------------------------- */
537
538 ASSERT (CHOLMOD(dump_sparse) (A, "A = triplet(T) result", Common) >= 0) ;
539 return (A) ;
540 }
541
542
543 /* ========================================================================== */
544 /* === cholmod_sparse_to_triplet ============================================ */
545 /* ========================================================================== */
546
547 /* Converts a sparse column-oriented matrix to triplet form.
548 * The resulting triplet matrix has the same xtype as the sparse matrix.
549 *
550 * workspace: none
551 */
552
CHOLMOD(sparse_to_triplet)553 cholmod_triplet *CHOLMOD(sparse_to_triplet)
554 (
555 /* ---- input ---- */
556 cholmod_sparse *A, /* matrix to copy */
557 /* --------------- */
558 cholmod_common *Common
559 )
560 {
561 double *Ax, *Az, *Tx, *Tz ;
562 Int *Ap, *Ai, *Ti, *Tj, *Anz ;
563 cholmod_triplet *T ;
564 Int i, xtype, p, pend, k, j, nrow, ncol, nz, stype, packed, up, lo,
565 both ;
566
567 /* ---------------------------------------------------------------------- */
568 /* check inputs */
569 /* ---------------------------------------------------------------------- */
570
571 RETURN_IF_NULL_COMMON (NULL) ;
572 RETURN_IF_NULL (A, NULL) ;
573 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
574 stype = SIGN (A->stype) ;
575 nrow = A->nrow ;
576 ncol = A->ncol ;
577 if (stype && nrow != ncol)
578 {
579 /* inputs invalid */
580 ERROR (CHOLMOD_INVALID, "matrix invalid") ;
581 return (NULL) ;
582 }
583 Ax = A->x ;
584 Az = A->z ;
585 xtype = A->xtype ;
586 Common->status = CHOLMOD_OK ;
587
588 ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;
589
590 /* ---------------------------------------------------------------------- */
591 /* allocate triplet matrix */
592 /* ---------------------------------------------------------------------- */
593
594 nz = CHOLMOD(nnz) (A, Common) ;
595 T = CHOLMOD(allocate_triplet) (nrow, ncol, nz, A->stype, A->xtype, Common) ;
596
597 if (Common->status < CHOLMOD_OK)
598 {
599 return (NULL) ; /* out of memory */
600 }
601
602 /* ---------------------------------------------------------------------- */
603 /* convert to a sparse matrix */
604 /* ---------------------------------------------------------------------- */
605
606 Ap = A->p ;
607 Ai = A->i ;
608 Anz = A->nz ;
609 packed = A->packed ;
610
611 Ti = T->i ;
612 Tj = T->j ;
613 Tx = T->x ;
614 Tz = T->z ;
615 T->stype = A->stype ;
616
617 both = (A->stype == 0) ;
618 up = (A->stype > 0) ;
619 lo = (A->stype < 0) ;
620
621 k = 0 ;
622
623 for (j = 0 ; j < ncol ; j++)
624 {
625 p = Ap [j] ;
626 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
627 for ( ; p < pend ; p++)
628 {
629 i = Ai [p] ;
630 if (both || (up && i <= j) || (lo && i >= j))
631 {
632 Ti [k] = Ai [p] ;
633 Tj [k] = j ;
634
635 if (xtype == CHOLMOD_REAL)
636 {
637 Tx [k] = Ax [p] ;
638 }
639 else if (xtype == CHOLMOD_COMPLEX)
640 {
641 Tx [2*k ] = Ax [2*p ] ;
642 Tx [2*k+1] = Ax [2*p+1] ;
643 }
644 else if (xtype == CHOLMOD_ZOMPLEX)
645 {
646 Tx [k] = Ax [p] ;
647 Tz [k] = Az [p] ;
648 }
649
650 k++ ;
651 ASSERT (k <= nz) ;
652 }
653 }
654 }
655
656 T->nnz = k ;
657
658 /* ---------------------------------------------------------------------- */
659 /* return result */
660 /* ---------------------------------------------------------------------- */
661
662 ASSERT (CHOLMOD(dump_triplet) (T, "T", Common)) ;
663 return (T) ;
664 }
665
666
667 /* ========================================================================== */
668 /* === cholmod_copy_triplet ================================================= */
669 /* ========================================================================== */
670
671 /* Create an exact copy of a triplet matrix, except that entries in unused
672 * space are not copied (they might not be initialized, and copying them would
673 * cause program checkers such as purify and valgrind to complain).
674 * The output triplet matrix has the same xtype as the input triplet matrix.
675 */
676
CHOLMOD(copy_triplet)677 cholmod_triplet *CHOLMOD(copy_triplet)
678 (
679 /* ---- input ---- */
680 cholmod_triplet *T, /* matrix to copy */
681 /* --------------- */
682 cholmod_common *Common
683 )
684 {
685 double *Tx, *Tz, *Cx, *Cz ;
686 Int *Ci, *Cj, *Ti, *Tj ;
687 cholmod_triplet *C ;
688 Int xtype, k, nz ;
689
690 /* ---------------------------------------------------------------------- */
691 /* check inputs */
692 /* ---------------------------------------------------------------------- */
693
694 RETURN_IF_NULL_COMMON (NULL) ;
695 RETURN_IF_NULL (T, NULL) ;
696 RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
697 nz = T->nnz ;
698 Ti = T->i ;
699 Tj = T->j ;
700 Tx = T->x ;
701 Tz = T->z ;
702 xtype = T->xtype ;
703 RETURN_IF_NULL (Ti, NULL) ;
704 RETURN_IF_NULL (Tj, NULL) ;
705 Common->status = CHOLMOD_OK ;
706 DEBUG (CHOLMOD(dump_triplet) (T, "T input", Common)) ;
707
708 /* ---------------------------------------------------------------------- */
709 /* allocate copy */
710 /* ---------------------------------------------------------------------- */
711
712 C = CHOLMOD(allocate_triplet) (T->nrow, T->ncol, T->nzmax, T->stype,
713 xtype, Common) ;
714
715 if (Common->status < CHOLMOD_OK)
716 {
717 return (NULL) ; /* out of memory */
718 }
719
720 /* ---------------------------------------------------------------------- */
721 /* copy the triplet matrix */
722 /* ---------------------------------------------------------------------- */
723
724 Ci = C->i ;
725 Cj = C->j ;
726 Cx = C->x ;
727 Cz = C->z ;
728 C->nnz = nz ;
729
730 for (k = 0 ; k < nz ; k++)
731 {
732 Ci [k] = Ti [k] ;
733 }
734 for (k = 0 ; k < nz ; k++)
735 {
736 Cj [k] = Tj [k] ;
737 }
738
739 if (xtype == CHOLMOD_REAL)
740 {
741 for (k = 0 ; k < nz ; k++)
742 {
743 Cx [k] = Tx [k] ;
744 }
745 }
746 else if (xtype == CHOLMOD_COMPLEX)
747 {
748 for (k = 0 ; k < nz ; k++)
749 {
750 Cx [2*k ] = Tx [2*k ] ;
751 Cx [2*k+1] = Tx [2*k+1] ;
752 }
753 }
754 else if (xtype == CHOLMOD_ZOMPLEX)
755 {
756 for (k = 0 ; k < nz ; k++)
757 {
758 Cx [k] = Tx [k] ;
759 Cz [k] = Tz [k] ;
760 }
761 }
762
763 /* ---------------------------------------------------------------------- */
764 /* return the result */
765 /* ---------------------------------------------------------------------- */
766
767 ASSERT (CHOLMOD(dump_triplet) (C, "C triplet copy", Common)) ;
768 return (C) ;
769 }
770