1 /* ========================================================================== */
2 /* === MatrixOps/cholmod_norm =============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/MatrixOps Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * http://www.suitesparse.com
8 * -------------------------------------------------------------------------- */
9
10 /* r = norm (A), compute the infinity-norm, 1-norm, or 2-norm of a sparse or
11 * dense matrix. Can compute the 2-norm only for a dense column vector.
12 * Returns -1 if an error occurs.
13 *
14 * Pattern, real, complex, and zomplex sparse matrices are supported.
15 */
16
17 #ifndef NGPL
18 #ifndef NMATRIXOPS
19
20 #include "cholmod_internal.h"
21 #include "cholmod_matrixops.h"
22
23
24 /* ========================================================================== */
25 /* === abs_value ============================================================ */
26 /* ========================================================================== */
27
28 /* Compute the absolute value of a real, complex, or zomplex value */
29
abs_value(int xtype,double * Ax,double * Az,Int p,cholmod_common * Common)30 static double abs_value
31 (
32 int xtype,
33 double *Ax,
34 double *Az,
35 Int p,
36 cholmod_common *Common
37 )
38 {
39 double s = 0 ;
40 switch (xtype)
41 {
42 case CHOLMOD_PATTERN:
43 s = 1 ;
44 break ;
45
46 case CHOLMOD_REAL:
47 s = fabs (Ax [p]) ;
48 break ;
49
50 case CHOLMOD_COMPLEX:
51 s = SuiteSparse_config.hypot_func (Ax [2*p], Ax [2*p+1]) ;
52 break ;
53
54 case CHOLMOD_ZOMPLEX:
55 s = SuiteSparse_config.hypot_func (Ax [p], Az [p]) ;
56 break ;
57 }
58 return (s) ;
59 }
60
61
62 /* ========================================================================== */
63 /* === cholmod_norm_dense =================================================== */
64 /* ========================================================================== */
65
CHOLMOD(norm_dense)66 double CHOLMOD(norm_dense)
67 (
68 /* ---- input ---- */
69 cholmod_dense *X, /* matrix to compute the norm of */
70 int norm, /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */
71 /* --------------- */
72 cholmod_common *Common
73 )
74 {
75 double xnorm, s, x, z ;
76 double *Xx, *Xz, *W ;
77 Int nrow, ncol, d, i, j, use_workspace, xtype ;
78
79 /* ---------------------------------------------------------------------- */
80 /* check inputs */
81 /* ---------------------------------------------------------------------- */
82
83 RETURN_IF_NULL_COMMON (EMPTY) ;
84 RETURN_IF_NULL (X, EMPTY) ;
85 RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, EMPTY) ;
86 Common->status = CHOLMOD_OK ;
87 ncol = X->ncol ;
88 if (norm < 0 || norm > 2 || (norm == 2 && ncol > 1))
89 {
90 ERROR (CHOLMOD_INVALID, "invalid norm") ;
91 return (EMPTY) ;
92 }
93
94 /* ---------------------------------------------------------------------- */
95 /* get inputs */
96 /* ---------------------------------------------------------------------- */
97
98 nrow = X->nrow ;
99 d = X->d ;
100 Xx = X->x ;
101 Xz = X->z ;
102 xtype = X->xtype ;
103
104 /* ---------------------------------------------------------------------- */
105 /* allocate workspace, if needed */
106 /* ---------------------------------------------------------------------- */
107
108 W = NULL ;
109 use_workspace = (norm == 0 && ncol > 4) ;
110 if (use_workspace)
111 {
112 CHOLMOD(allocate_work) (0, 0, nrow, Common) ;
113 W = Common->Xwork ;
114 if (Common->status < CHOLMOD_OK)
115 {
116 /* oops, no workspace */
117 use_workspace = FALSE ;
118 }
119 }
120
121
122 /* ---------------------------------------------------------------------- */
123 /* compute the norm */
124 /* ---------------------------------------------------------------------- */
125
126 xnorm = 0 ;
127
128 if (use_workspace)
129 {
130
131 /* ------------------------------------------------------------------ */
132 /* infinity-norm = max row sum, using stride-1 access of X */
133 /* ------------------------------------------------------------------ */
134
135 DEBUG (for (i = 0 ; i < nrow ; i++) ASSERT (W [i] == 0)) ;
136
137 /* this is faster than stride-d, but requires O(nrow) workspace */
138 for (j = 0 ; j < ncol ; j++)
139 {
140 for (i = 0 ; i < nrow ; i++)
141 {
142 W [i] += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
143 }
144 }
145 for (i = 0 ; i < nrow ; i++)
146 {
147 s = W [i] ;
148 if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
149 {
150 xnorm = s ;
151 }
152 W [i] = 0 ;
153 }
154
155 }
156 else if (norm == 0)
157 {
158
159 /* ------------------------------------------------------------------ */
160 /* infinity-norm = max row sum, using stride-d access of X */
161 /* ------------------------------------------------------------------ */
162
163 for (i = 0 ; i < nrow ; i++)
164 {
165 s = 0 ;
166 for (j = 0 ; j < ncol ; j++)
167 {
168 s += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
169 }
170 if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
171 {
172 xnorm = s ;
173 }
174 }
175
176 }
177 else if (norm == 1)
178 {
179
180 /* ------------------------------------------------------------------ */
181 /* 1-norm = max column sum */
182 /* ------------------------------------------------------------------ */
183
184 for (j = 0 ; j < ncol ; j++)
185 {
186 s = 0 ;
187 for (i = 0 ; i < nrow ; i++)
188 {
189 s += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
190 }
191 if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
192 {
193 xnorm = s ;
194 }
195 }
196 }
197 else
198 {
199
200 /* ------------------------------------------------------------------ */
201 /* 2-norm = sqrt (sum (X.^2)) */
202 /* ------------------------------------------------------------------ */
203
204 switch (xtype)
205 {
206
207 case CHOLMOD_REAL:
208 for (i = 0 ; i < nrow ; i++)
209 {
210 x = Xx [i] ;
211 xnorm += x*x ;
212 }
213 break ;
214
215 case CHOLMOD_COMPLEX:
216 for (i = 0 ; i < nrow ; i++)
217 {
218 x = Xx [2*i ] ;
219 z = Xx [2*i+1] ;
220 xnorm += x*x + z*z ;
221 }
222 break ;
223
224 case CHOLMOD_ZOMPLEX:
225 for (i = 0 ; i < nrow ; i++)
226 {
227 x = Xx [i] ;
228 z = Xz [i] ;
229 xnorm += x*x + z*z ;
230 }
231 break ;
232 }
233
234 xnorm = sqrt (xnorm) ;
235 }
236
237 /* ---------------------------------------------------------------------- */
238 /* return result */
239 /* ---------------------------------------------------------------------- */
240
241 return (xnorm) ;
242 }
243
244
245 /* ========================================================================== */
246 /* === cholmod_norm_sparse ================================================== */
247 /* ========================================================================== */
248
CHOLMOD(norm_sparse)249 double CHOLMOD(norm_sparse)
250 (
251 /* ---- input ---- */
252 cholmod_sparse *A, /* matrix to compute the norm of */
253 int norm, /* type of norm: 0: inf. norm, 1: 1-norm */
254 /* --------------- */
255 cholmod_common *Common
256 )
257 {
258 double anorm, s ;
259 double *Ax, *Az, *W ;
260 Int *Ap, *Ai, *Anz ;
261 Int i, j, p, pend, nrow, ncol, packed, xtype ;
262
263 /* ---------------------------------------------------------------------- */
264 /* check inputs */
265 /* ---------------------------------------------------------------------- */
266
267 RETURN_IF_NULL_COMMON (EMPTY) ;
268 RETURN_IF_NULL (A, EMPTY) ;
269 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
270 Common->status = CHOLMOD_OK ;
271 ncol = A->ncol ;
272 nrow = A->nrow ;
273 if (norm < 0 || norm > 1)
274 {
275 ERROR (CHOLMOD_INVALID, "invalid norm") ;
276 return (EMPTY) ;
277 }
278 if (A->stype && nrow != ncol)
279 {
280 ERROR (CHOLMOD_INVALID, "matrix invalid") ;
281 return (EMPTY) ;
282 }
283
284 /* ---------------------------------------------------------------------- */
285 /* get inputs */
286 /* ---------------------------------------------------------------------- */
287
288 Ap = A->p ;
289 Ai = A->i ;
290 Ax = A->x ;
291 Az = A->z ;
292 Anz = A->nz ;
293 packed = A->packed ;
294 xtype = A->xtype ;
295
296 /* ---------------------------------------------------------------------- */
297 /* allocate workspace, if needed */
298 /* ---------------------------------------------------------------------- */
299
300 W = NULL ;
301 if (A->stype || norm == 0)
302 {
303 CHOLMOD(allocate_work) (0, 0, nrow, Common) ;
304 W = Common->Xwork ;
305 if (Common->status < CHOLMOD_OK)
306 {
307 /* out of memory */
308 return (EMPTY) ;
309 }
310 DEBUG (for (i = 0 ; i < nrow ; i++) ASSERT (W [i] == 0)) ;
311 }
312
313 /* ---------------------------------------------------------------------- */
314 /* compute the norm */
315 /* ---------------------------------------------------------------------- */
316
317 anorm = 0 ;
318
319 if (A->stype > 0)
320 {
321
322 /* ------------------------------------------------------------------ */
323 /* A is symmetric with upper triangular part stored */
324 /* ------------------------------------------------------------------ */
325
326 /* infinity-norm = 1-norm = max row/col sum */
327 for (j = 0 ; j < ncol ; j++)
328 {
329 p = Ap [j] ;
330 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
331 for ( ; p < pend ; p++)
332 {
333 i = Ai [p] ;
334 s = abs_value (xtype, Ax, Az, p, Common) ;
335 if (i == j)
336 {
337 W [i] += s ;
338 }
339 else if (i < j)
340 {
341 W [i] += s ;
342 W [j] += s ;
343 }
344 }
345 }
346
347 }
348 else if (A->stype < 0)
349 {
350
351 /* ------------------------------------------------------------------ */
352 /* A is symmetric with lower triangular part stored */
353 /* ------------------------------------------------------------------ */
354
355 /* infinity-norm = 1-norm = max row/col sum */
356 for (j = 0 ; j < ncol ; j++)
357 {
358 p = Ap [j] ;
359 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
360 for ( ; p < pend ; p++)
361 {
362 i = Ai [p] ;
363 s = abs_value (xtype, Ax, Az, p, Common) ;
364 if (i == j)
365 {
366 W [i] += s ;
367 }
368 else if (i > j)
369 {
370 W [i] += s ;
371 W [j] += s ;
372 }
373 }
374 }
375
376 }
377 else if (norm == 0)
378 {
379
380 /* ------------------------------------------------------------------ */
381 /* A is unsymmetric, compute the infinity-norm */
382 /* ------------------------------------------------------------------ */
383
384 /* infinity-norm = max row sum */
385 for (j = 0 ; j < ncol ; j++)
386 {
387 p = Ap [j] ;
388 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
389 for ( ; p < pend ; p++)
390 {
391 W [Ai [p]] += abs_value (xtype, Ax, Az, p, Common) ;
392 }
393 }
394
395 }
396 else
397 {
398
399 /* ------------------------------------------------------------------ */
400 /* A is unsymmetric, compute the 1-norm */
401 /* ------------------------------------------------------------------ */
402
403 /* 1-norm = max column sum */
404 for (j = 0 ; j < ncol ; j++)
405 {
406 p = Ap [j] ;
407 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
408 if (xtype == CHOLMOD_PATTERN)
409 {
410 s = pend - p ;
411 }
412 else
413 {
414 s = 0 ;
415 for ( ; p < pend ; p++)
416 {
417 s += abs_value (xtype, Ax, Az, p, Common) ;
418 }
419 }
420 if ((IS_NAN (s) || s > anorm) && !IS_NAN (anorm))
421 {
422 anorm = s ;
423 }
424 }
425 }
426
427 /* ---------------------------------------------------------------------- */
428 /* compute the max row sum */
429 /* ---------------------------------------------------------------------- */
430
431 if (A->stype || norm == 0)
432 {
433 for (i = 0 ; i < nrow ; i++)
434 {
435 s = W [i] ;
436 if ((IS_NAN (s) || s > anorm) && !IS_NAN (anorm))
437 {
438 anorm = s ;
439 }
440 W [i] = 0 ;
441 }
442 }
443
444 /* ---------------------------------------------------------------------- */
445 /* return result */
446 /* ---------------------------------------------------------------------- */
447
448 return (anorm) ;
449 }
450 #endif
451 #endif
452