1 /*! \file
2 Copyright (c) 2003, The Regents of the University of California, through
3 Lawrence Berkeley National Laboratory (subject to receipt of any required
4 approvals from U.S. Dept. of Energy)
5
6 All rights reserved.
7
8 The source code is distributed under BSD license, see the file License.txt
9 at the top-level directory.
10 */
11
12 /*! @file cdiagonal.c
13 * \brief Auxiliary routines to work with diagonal elements
14 *
15 * <pre>
16 * -- SuperLU routine (version 4.0) --
17 * Lawrence Berkeley National Laboratory
18 * June 30, 2009
19 * </pre>
20 */
21
22 #include "slu_cdefs.h"
23
cfill_diag(int n,NCformat * Astore)24 int cfill_diag(int n, NCformat *Astore)
25 /* fill explicit zeros on the diagonal entries, so that the matrix is not
26 structurally singular. */
27 {
28 complex *nzval = (complex *)Astore->nzval;
29 int *rowind = Astore->rowind;
30 int *colptr = Astore->colptr;
31 int nnz = colptr[n];
32 int fill = 0;
33 complex *nzval_new;
34 complex zero = {0.0, 0.0};
35 int *rowind_new;
36 int i, j, diag;
37
38 for (i = 0; i < n; i++)
39 {
40 diag = -1;
41 for (j = colptr[i]; j < colptr[i + 1]; j++)
42 if (rowind[j] == i) diag = j;
43 if (diag < 0) fill++;
44 }
45 if (fill)
46 {
47 nzval_new = complexMalloc(nnz + fill);
48 rowind_new = intMalloc(nnz + fill);
49 fill = 0;
50 for (i = 0; i < n; i++)
51 {
52 diag = -1;
53 for (j = colptr[i] - fill; j < colptr[i + 1]; j++)
54 {
55 if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;
56 nzval_new[j + fill] = nzval[j];
57 }
58 if (diag < 0)
59 {
60 rowind_new[colptr[i + 1] + fill] = i;
61 nzval_new[colptr[i + 1] + fill] = zero;
62 fill++;
63 }
64 colptr[i + 1] += fill;
65 }
66 Astore->nzval = nzval_new;
67 Astore->rowind = rowind_new;
68 SUPERLU_FREE(nzval);
69 SUPERLU_FREE(rowind);
70 }
71 Astore->nnz += fill;
72 return fill;
73 }
74
cdominate(int n,NCformat * Astore)75 int cdominate(int n, NCformat *Astore)
76 /* make the matrix diagonally dominant */
77 {
78 complex *nzval = (complex *)Astore->nzval;
79 int *rowind = Astore->rowind;
80 int *colptr = Astore->colptr;
81 int nnz = colptr[n];
82 int fill = 0;
83 complex *nzval_new;
84 int *rowind_new;
85 int i, j, diag;
86 double s;
87
88 for (i = 0; i < n; i++)
89 {
90 diag = -1;
91 for (j = colptr[i]; j < colptr[i + 1]; j++)
92 if (rowind[j] == i) diag = j;
93 if (diag < 0) fill++;
94 }
95 if (fill)
96 {
97 nzval_new = complexMalloc(nnz + fill);
98 rowind_new = intMalloc(nnz+ fill);
99 fill = 0;
100 for (i = 0; i < n; i++)
101 {
102 s = 1e-6;
103 diag = -1;
104 for (j = colptr[i] - fill; j < colptr[i + 1]; j++)
105 {
106 if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;
107 nzval_new[j + fill] = nzval[j];
108 s += c_abs1(&nzval_new[j + fill]);
109 }
110 if (diag >= 0) {
111 nzval_new[diag+fill].r = s * 3.0;
112 nzval_new[diag+fill].i = 0.0;
113 } else {
114 rowind_new[colptr[i + 1] + fill] = i;
115 nzval_new[colptr[i + 1] + fill].r = s * 3.0;
116 nzval_new[colptr[i + 1] + fill].i = 0.0;
117 fill++;
118 }
119 colptr[i + 1] += fill;
120 }
121 Astore->nzval = nzval_new;
122 Astore->rowind = rowind_new;
123 SUPERLU_FREE(nzval);
124 SUPERLU_FREE(rowind);
125 }
126 else
127 {
128 for (i = 0; i < n; i++)
129 {
130 s = 1e-6;
131 diag = -1;
132 for (j = colptr[i]; j < colptr[i + 1]; j++)
133 {
134 if (rowind[j] == i) diag = j;
135 s += c_abs1(&nzval[j]);
136 }
137 nzval[diag].r = s * 3.0;
138 nzval[diag].i = 0.0;
139 }
140 }
141 Astore->nnz += fill;
142 return fill;
143 }
144