Amesos2 - Direct Sparse Solver Interfaces Version of the Day
klu2_scale.hpp
1/* ========================================================================== */
2/* === KLU_scale ============================================================ */
3/* ========================================================================== */
4// @HEADER
5// ***********************************************************************
6//
7// KLU2: A Direct Linear Solver package
8// Copyright 2011 Sandia Corporation
9//
10// Under terms of Contract DE-AC04-94AL85000, with Sandia Corporation, the
11// U.S. Government retains certain rights in this software.
12//
13// This library is free software; you can redistribute it and/or modify
14// it under the terms of the GNU Lesser General Public License as
15// published by the Free Software Foundation; either version 2.1 of the
16// License, or (at your option) any later version.
17//
18// This library is distributed in the hope that it will be useful, but
19// WITHOUT ANY WARRANTY; without even the implied warranty of
20// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21// Lesser General Public License for more details.
22//
23// You should have received a copy of the GNU Lesser General Public
24// License along with this library; if not, write to the Free Software
25// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
26// USA
27// Questions? Contact Mike A. Heroux (maherou@sandia.gov)
28//
29// KLU2 is derived work from KLU, licensed under LGPL, and copyrighted by
30// University of Florida. The Authors of KLU are Timothy A. Davis and
31// Eka Palamadai. See Doc/KLU_README.txt for the licensing and copyright
32// information for KLU.
33//
34// ***********************************************************************
35// @HEADER
36
37/* Scale a matrix and check to see if it is valid. Can be called by the user.
38 * This is called by KLU_factor and KLU_refactor. Returns TRUE if the input
39 * matrix is valid, FALSE otherwise. If the W input argument is non-NULL,
40 * then the input matrix is checked for duplicate entries.
41 *
42 * scaling methods:
43 * <0: no scaling, do not compute Rs, and do not check input matrix.
44 * 0: no scaling
45 * 1: the scale factor for row i is sum (abs (A (i,:)))
46 * 2 or more: the scale factor for row i is max (abs (A (i,:)))
47 */
48
49#ifndef KLU2_SCALE_HPP
50#define KLU2_SCALE_HPP
51
52#include "klu2_internal.h"
53
54template <typename Entry, typename Int>
55Int KLU_scale /* return TRUE if successful, FALSE otherwise */
56(
57 /* inputs, not modified */
58 Int scale, /* 0: none, 1: sum, 2: max */
59 Int n,
60 Int Ap [ ], /* size n+1, column pointers */
61 Int Ai [ ], /* size nz, row indices */
62 /* TODO : double to Entry */
63 Entry Ax [ ],
64 /* outputs, not defined on input */
65 /* TODO : double to Entry */
66 double Rs [ ], /* size n, can be NULL if scale <= 0 */
67 /* workspace, not defined on input or output */
68 Int W [ ], /* size n, can be NULL */
69 /* --------------- */
70 KLU_common<Entry, Int> *Common
71)
72{
73 double a ;
74 Entry *Az ;
75 Int row, col, p, pend, check_duplicates ;
76
77 /* ---------------------------------------------------------------------- */
78 /* check inputs */
79 /* ---------------------------------------------------------------------- */
80
81 if (Common == NULL)
82 {
83 return (FALSE) ;
84 }
85 Common->status = KLU_OK ;
86
87 if (scale < 0)
88 {
89 /* return without checking anything and without computing the
90 * scale factors */
91 return (TRUE) ;
92 }
93
94 Az = (Entry *) Ax ;
95
96 if (n <= 0 || Ap == NULL || Ai == NULL || Az == NULL ||
97 (scale > 0 && Rs == NULL))
98 {
99 /* Ap, Ai, Ax and Rs must be present, and n must be > 0 */
100 Common->status = KLU_INVALID ;
101 return (FALSE) ;
102 }
103 if (Ap [0] != 0 || Ap [n] < 0)
104 {
105 /* nz = Ap [n] must be >= 0 and Ap [0] must equal zero */
106 Common->status = KLU_INVALID ;
107 return (FALSE) ;
108 }
109 for (col = 0 ; col < n ; col++)
110 {
111 if (Ap [col] > Ap [col+1])
112 {
113 /* column pointers must be non-decreasing */
114 Common->status = KLU_INVALID ;
115 return (FALSE) ;
116 }
117 }
118
119 /* ---------------------------------------------------------------------- */
120 /* scale */
121 /* ---------------------------------------------------------------------- */
122
123 if (scale > 0)
124 {
125 /* initialize row sum or row max */
126 for (row = 0 ; row < n ; row++)
127 {
128 Rs [row] = 0 ;
129 }
130 }
131
132 /* check for duplicates only if W is present */
133 check_duplicates = (W != (Int *) NULL) ;
134 if (check_duplicates)
135 {
136 for (row = 0 ; row < n ; row++)
137 {
138 W [row] = AMESOS2_KLU2_EMPTY ;
139 }
140 }
141
142 for (col = 0 ; col < n ; col++)
143 {
144 pend = Ap [col+1] ;
145 for (p = Ap [col] ; p < pend ; p++)
146 {
147 row = Ai [p] ;
148 if (row < 0 || row >= n)
149 {
150 /* row index out of range, or duplicate entry */
151 Common->status = KLU_INVALID ;
152 return (FALSE) ;
153 }
154 if (check_duplicates)
155 {
156 if (W [row] == col)
157 {
158 /* duplicate entry */
159 Common->status = KLU_INVALID ;
160 return (FALSE) ;
161 }
162 /* flag row i as appearing in column col */
163 W [row] = col ;
164 }
165 /* a = ABS (Az [p]) ;*/
166 KLU2_ABS (a, Az [p]) ;
167 if (scale == 1)
168 {
169 /* accumulate the abs. row sum */
170 Rs [row] += a ;
171 }
172 else if (scale > 1)
173 {
174 /* find the max abs. value in the row */
175 Rs [row] = MAX (Rs [row], a) ;
176 }
177 }
178 }
179
180 if (scale > 0)
181 {
182 /* do not scale empty rows */
183 for (row = 0 ; row < n ; row++)
184 {
185 /* matrix is singular */
186 PRINTF (("Rs [%d] = %g\n", row, Rs [row])) ;
187
188 if (Rs [row] == 0.0)
189 {
190 PRINTF (("Row %d of A is all zero\n", row)) ;
191 Rs [row] = 1.0 ;
192 }
193 }
194 }
195
196 return (TRUE) ;
197}
198#endif
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.