← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to scipy/sparse/linalg/dsolve/SuperLU/SRC/cdiagonal.c

  • Committer: Bazaar Package Importer
  • Author(s): Sameer Morar
  • Date: 2011-02-03 04:28:09 UTC
  • mfrom: (9.1.2 experimental)
  • Revision ID: james.westby@ubuntu.com-20110203042809-qs95kuod7723re62
Tags: 0.8.0+dfsg1-1ubuntu1
https://launchpad.net/bugs/696403
* Merge from debian experimental (LP: #696403). Remaining changes:
  - debian/patches/stdc_format_macros.patch: Fix FTBFS issue with python 2.7

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/sparse/linalg/dsolve/SuperLU/SRC/cdiagonal.c
 
1
 
 
2
/*! @file cdiagonal.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_cdefs.h"
 
13
 
 
14
int cfill_diag(int n, NCformat *Astore)
 
15
/* fill explicit zeros on the diagonal entries, so that the matrix is not
 
16
   structurally singular. */
 
17
{
 
18
    complex *nzval = (complex *)Astore->nzval;
 
19
    int *rowind = Astore->rowind;
 
20
    int *colptr = Astore->colptr;
 
21
    int nnz = colptr[n];
 
22
    int fill = 0;
 
23
    complex *nzval_new;
 
24
    complex zero = {1.0, 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 = complexMalloc(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
 
 
65
int cdominate(int n, NCformat *Astore)
 
66
/* make the matrix diagonally dominant */
 
67
{
 
68
    complex *nzval = (complex *)Astore->nzval;
 
69
    int *rowind = Astore->rowind;
 
70
    int *colptr = Astore->colptr;
 
71
    int nnz = colptr[n];
 
72
    int fill = 0;
 
73
    complex *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 = complexMalloc(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
                nzval_new[j + fill] = nzval[j];
 
98
                s += slu_c_abs1(&nzval_new[j + fill]);
 
99
            }
 
100
            if (diag >= 0) {
 
101
                nzval_new[diag+fill].r = s * 3.0;
 
102
                nzval_new[diag+fill].i = 0.0;
 
103
            } else {
 
104
                rowind_new[colptr[i + 1] + fill] = i;
 
105
                nzval_new[colptr[i + 1] + fill].r = s * 3.0;
 
106
                nzval_new[colptr[i + 1] + fill].i = 0.0;
 
107
                fill++;
 
108
            }
 
109
            colptr[i + 1] += fill;
 
110
        }
 
111
        Astore->nzval = nzval_new;
 
112
        Astore->rowind = rowind_new;
 
113
        SUPERLU_FREE(nzval);
 
114
        SUPERLU_FREE(rowind);
 
115
    }
 
116
    else
 
117
    {
 
118
        for (i = 0; i < n; i++)
 
119
        {
 
120
            s = 1e-6;
 
121
            diag = -1;
 
122
            for (j = colptr[i]; j < colptr[i + 1]; j++)
 
123
            {
 
124
                if (rowind[j] == i) diag = j;
 
125
                s += slu_c_abs1(&nzval[j]);
 
126
            }
 
127
            nzval[diag].r = s * 3.0;
 
128
            nzval[diag].i = 0.0;
 
129
        }
 
130
    }
 
131
    Astore->nnz += fill;
 
132
    return fill;
 
133
}