1 /* ========================================================================== */
2 /* === MatrixOps/cholmod_scale ============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MatrixOps Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* scale a matrix: A = diag(s)*A, A*diag(s), s*A, or diag(s)*A*diag(s)
11 *
12 * A can be of any type (packed/unpacked, upper/lower/unsymmetric).
13 * The symmetry of A is ignored; all entries in the matrix are modified.
14 *
15 * If A is m-by-n unsymmetric but scaled symmtrically, the result is
16 * A = diag (s (1:m)) * A * diag (s (1:n)).
17 *
18 * Note: diag(s) should be interpretted as spdiags(s,0,n,n) where n=length(s).
19 *
20 * Row or column scaling of a symmetric matrix still results in a symmetric
21 * matrix, since entries are still ignored by other routines.
22 * For example, when row-scaling a symmetric matrix where just the upper
23 * triangular part is stored (and lower triangular entries ignored)
24 * A = diag(s)*triu(A) is performed, where the result A is also
25 * symmetric-upper. This has the effect of modifying the implicit lower
26 * triangular part. In MATLAB notation:
27 *
28 * U = diag(s)*triu(A) ;
29 * L = tril (U',-1)
30 * A = L + U ;
31 *
32 * The scale parameter determines the kind of scaling to perform:
33 *
34 * CHOLMOD_SCALAR: s[0]*A
35 * CHOLMOD_ROW: diag(s)*A
36 * CHOLMOD_COL: A*diag(s)
37 * CHOLMOD_SYM: diag(s)*A*diag(s)
38 *
39 * The size of S depends on the scale parameter:
40 *
41 * CHOLMOD_SCALAR: size 1
42 * CHOLMOD_ROW: size nrow-by-1 or 1-by-nrow
43 * CHOLMOD_COL: size ncol-by-1 or 1-by-ncol
44 * CHOLMOD_SYM: size max(nrow,ncol)-by-1, or 1-by-max(nrow,ncol)
45 *
46 * workspace: none
47 *
48 * Only real matrices are supported.
49 */
50
51 #ifndef NGPL
52 #ifndef NMATRIXOPS
53
54 #include "cholmod_internal.h"
55 #include "cholmod_matrixops.h"
56
57
58 /* ========================================================================== */
59 /* === cholmod_scale ======================================================== */
60 /* ========================================================================== */
61
CHOLMOD(scale)62 int CHOLMOD(scale)
63 (
64 /* ---- input ---- */
65 cholmod_dense *S, /* scale factors (scalar or vector) */
66 int scale, /* type of scaling to compute */
67 /* ---- in/out --- */
68 cholmod_sparse *A, /* matrix to scale */
69 /* --------------- */
70 cholmod_common *Common
71 )
72 {
73 double t ;
74 double *Ax, *s ;
75 Int *Ap, *Anz, *Ai ;
76 Int packed, j, ncol, nrow, p, pend, sncol, snrow, nn, ok ;
77
78 /* ---------------------------------------------------------------------- */
79 /* check inputs */
80 /* ---------------------------------------------------------------------- */
81
82 RETURN_IF_NULL_COMMON (FALSE) ;
83 RETURN_IF_NULL (A, FALSE) ;
84 RETURN_IF_NULL (S, FALSE) ;
85 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
86 RETURN_IF_XTYPE_INVALID (S, CHOLMOD_REAL, CHOLMOD_REAL, FALSE) ;
87 ncol = A->ncol ;
88 nrow = A->nrow ;
89 sncol = S->ncol ;
90 snrow = S->nrow ;
91 if (scale == CHOLMOD_SCALAR)
92 {
93 ok = (snrow == 1 && sncol == 1) ;
94 }
95 else if (scale == CHOLMOD_ROW)
96 {
97 ok = (snrow == nrow && sncol == 1) || (snrow == 1 && sncol == nrow) ;
98 }
99 else if (scale == CHOLMOD_COL)
100 {
101 ok = (snrow == ncol && sncol == 1) || (snrow == 1 && sncol == ncol) ;
102 }
103 else if (scale == CHOLMOD_SYM)
104 {
105 nn = MAX (nrow, ncol) ;
106 ok = (snrow == nn && sncol == 1) || (snrow == 1 && sncol == nn) ;
107 }
108 else
109 {
110 /* scale invalid */
111 ERROR (CHOLMOD_INVALID, "invalid scaling option") ;
112 return (FALSE) ;
113 }
114 if (!ok)
115 {
116 /* S is wrong size */
117 ERROR (CHOLMOD_INVALID, "invalid scale factors") ;
118 return (FALSE) ;
119 }
120 Common->status = CHOLMOD_OK ;
121
122 /* ---------------------------------------------------------------------- */
123 /* get inputs */
124 /* ---------------------------------------------------------------------- */
125
126 Ap = A->p ;
127 Anz = A->nz ;
128 Ai = A->i ;
129 Ax = A->x ;
130 packed = A->packed ;
131 s = S->x ;
132
133 /* ---------------------------------------------------------------------- */
134 /* scale the matrix */
135 /* ---------------------------------------------------------------------- */
136
137 if (scale == CHOLMOD_ROW)
138 {
139
140 /* ------------------------------------------------------------------ */
141 /* A = diag(s)*A, row scaling */
142 /* ------------------------------------------------------------------ */
143
144 for (j = 0 ; j < ncol ; j++)
145 {
146 p = Ap [j] ;
147 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
148 for ( ; p < pend ; p++)
149 {
150 Ax [p] *= s [Ai [p]] ;
151 }
152 }
153
154 }
155 else if (scale == CHOLMOD_COL)
156 {
157
158 /* ------------------------------------------------------------------ */
159 /* A = A*diag(s), column scaling */
160 /* ------------------------------------------------------------------ */
161
162 for (j = 0 ; j < ncol ; j++)
163 {
164 t = s [j] ;
165 p = Ap [j] ;
166 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
167 for ( ; p < pend ; p++)
168 {
169 Ax [p] *= t ;
170 }
171 }
172
173 }
174 else if (scale == CHOLMOD_SYM)
175 {
176
177 /* ------------------------------------------------------------------ */
178 /* A = diag(s)*A*diag(s), symmetric scaling */
179 /* ------------------------------------------------------------------ */
180
181 for (j = 0 ; j < ncol ; j++)
182 {
183 t = s [j] ;
184 p = Ap [j] ;
185 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
186 for ( ; p < pend ; p++)
187 {
188 Ax [p] *= t * s [Ai [p]] ;
189 }
190 }
191
192 }
193 else if (scale == CHOLMOD_SCALAR)
194 {
195
196 /* ------------------------------------------------------------------ */
197 /* A = s[0] * A, scalar scaling */
198 /* ------------------------------------------------------------------ */
199
200 t = s [0] ;
201 for (j = 0 ; j < ncol ; j++)
202 {
203 p = Ap [j] ;
204 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
205 for ( ; p < pend ; p++)
206 {
207 Ax [p] *= t ;
208 }
209 }
210 }
211
212 ASSERT (CHOLMOD(dump_sparse) (A, "A scaled", Common) >= 0) ;
213 return (TRUE) ;
214 }
215 #endif
216 #endif
217