1
2 /*! @file sdiagonal.c
3 * \brief Auxiliary routines to work with diagonal elements
4 *
5 * <pre>
6 * -- SuperLU routine (version 4.0) --
7 * Lawrence Berkeley National Laboratory
8 * June 30, 2009
9 * </pre>
10 */
11
12 #include "slu_sdefs.h"
13
sfill_diag(int n,NCformat * Astore)14 int sfill_diag(int n, NCformat *Astore)
15 /* fill explicit zeros on the diagonal entries, so that the matrix is not
16 structurally singular. */
17 {
18 float *nzval = (float *)Astore->nzval;
19 int *rowind = Astore->rowind;
20 int *colptr = Astore->colptr;
21 int nnz = colptr[n];
22 int fill = 0;
23 float *nzval_new;
24 float zero = 0.0;
25 int *rowind_new;
26 int i, j, diag;
27
28 for (i = 0; i < n; i++)
29 {
30 diag = -1;
31 for (j = colptr[i]; j < colptr[i + 1]; j++)
32 if (rowind[j] == i) diag = j;
33 if (diag < 0) fill++;
34 }
35 if (fill)
36 {
37 nzval_new = floatMalloc(nnz + fill);
38 rowind_new = intMalloc(nnz + fill);
39 fill = 0;
40 for (i = 0; i < n; i++)
41 {
42 diag = -1;
43 for (j = colptr[i] - fill; j < colptr[i + 1]; j++)
44 {
45 if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;
46 nzval_new[j + fill] = nzval[j];
47 }
48 if (diag < 0)
49 {
50 rowind_new[colptr[i + 1] + fill] = i;
51 nzval_new[colptr[i + 1] + fill] = zero;
52 fill++;
53 }
54 colptr[i + 1] += fill;
55 }
56 Astore->nzval = nzval_new;
57 Astore->rowind = rowind_new;
58 SUPERLU_FREE(nzval);
59 SUPERLU_FREE(rowind);
60 }
61 Astore->nnz += fill;
62 return fill;
63 }
64
sdominate(int n,NCformat * Astore)65 int sdominate(int n, NCformat *Astore)
66 /* make the matrix diagonally dominant */
67 {
68 float *nzval = (float *)Astore->nzval;
69 int *rowind = Astore->rowind;
70 int *colptr = Astore->colptr;
71 int nnz = colptr[n];
72 int fill = 0;
73 float *nzval_new;
74 int *rowind_new;
75 int i, j, diag;
76 double s;
77
78 for (i = 0; i < n; i++)
79 {
80 diag = -1;
81 for (j = colptr[i]; j < colptr[i + 1]; j++)
82 if (rowind[j] == i) diag = j;
83 if (diag < 0) fill++;
84 }
85 if (fill)
86 {
87 nzval_new = floatMalloc(nnz + fill);
88 rowind_new = intMalloc(nnz+ fill);
89 fill = 0;
90 for (i = 0; i < n; i++)
91 {
92 s = 1e-6;
93 diag = -1;
94 for (j = colptr[i] - fill; j < colptr[i + 1]; j++)
95 {
96 if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;
97 s += fabs(nzval_new[j + fill] = nzval[j]);
98 }
99 if (diag >= 0) {
100 nzval_new[diag+fill] = s * 3.0;
101 } else {
102 rowind_new[colptr[i + 1] + fill] = i;
103 nzval_new[colptr[i + 1] + fill] = s * 3.0;
104 fill++;
105 }
106 colptr[i + 1] += fill;
107 }
108 Astore->nzval = nzval_new;
109 Astore->rowind = rowind_new;
110 SUPERLU_FREE(nzval);
111 SUPERLU_FREE(rowind);
112 }
113 else
114 {
115 for (i = 0; i < n; i++)
116 {
117 s = 1e-6;
118 diag = -1;
119 for (j = colptr[i]; j < colptr[i + 1]; j++)
120 {
121 if (rowind[j] == i) diag = j;
122 s += fabs(nzval[j]);
123 }
124 nzval[diag] = s * 3.0;
125 }
126 }
127 Astore->nnz += fill;
128 return fill;
129 }
130