basiclu-sys 0.1.2

Rust FFI bindings for BasicLU
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

/*
 * lu_solve_triangular.c
 *
 * Copyright (C) 2016-2018  ERGO-Code
 *
 */

#include "lu_internal.h"


/**
 * lu_solve_triangular() - substitution with triangular matrix
 *
 * The symbolic nonzero pattern of the solution must be given in topological
 * order in pattern_symb[0..nz_symb-1]. On return pattern[0..nz-1] holds the
 * nonzero pattern of the solution after dropping numerical zeros; nz is
 * returned. pattern and pattern_symb can point to the same array.
 *
 * Entries in the solution that are less than or equal to droptol are set to
 * zero. When droptol is zero or negative, then no entries will be set to zero.
 *
 * Note: The nonzero pattern of the solution never includes zeros. That means,
 *       even if droptol is negative, the output pattern is not identical to
 *       the symbolic pattern when exact cancellation happens.
 *
 * The pivot elements must be stored separately to the matrix. When pivot is
 * NULL, then the pivot elements are assumed to be 1. The matrix is given in
 * parallel arrays index, value. When end is not NULL, column j has elements
 *
 *   index[begin[j]..end[j]-1], value[begin[j]..end[j]-1].
 *
 * When end is NULL, then each column must be terminated by a negative index.
 *
 */
lu_int lu_solve_triangular
(
    const lu_int nz_symb,
    const lu_int *pattern_symb,
    const lu_int *begin,
    const lu_int *end,
    const lu_int *index,
    const double *value,
    const double *pivot,
    const double droptol,
    double *lhs,                /* solution overwrites RHS */
    lu_int *pattern,
    lu_int *flops               /* add flop count */
)
{
    lu_int i, ipivot, pos, n, nz = 0, flop_count = 0;
    double x;

    if (pivot && end)
    {
        for (n = 0; n < nz_symb; n++)
        {
            ipivot = pattern_symb[n];
            if (lhs[ipivot])
            {
                x = lhs[ipivot] /= pivot[ipivot];
                flop_count++;
                for (pos = begin[ipivot]; pos < end[ipivot]; pos++)
                {
                    i = index[pos];
                    lhs[i] -= x * value[pos];
                    flop_count++;
                }
                if (fabs(x) > droptol)
                    pattern[nz++] = ipivot;
                else
                    lhs[ipivot] = 0.0;
            }
        }
    }
    else if (pivot)
    {
        for (n = 0; n < nz_symb; n++)
        {
            ipivot = pattern_symb[n];
            if (lhs[ipivot])
            {
                x = lhs[ipivot] /= pivot[ipivot];
                flop_count++;
                for (pos = begin[ipivot]; (i = index[pos]) >= 0; pos++)
                {
                    lhs[i] -= x * value[pos];
                    flop_count++;
                }
                if (fabs(x) > droptol)
                    pattern[nz++] = ipivot;
                else
                    lhs[ipivot] = 0.0;
            }
        }
    }
    else if (end)
    {
        for (n = 0; n < nz_symb; n++)
        {
            ipivot = pattern_symb[n];
            if (lhs[ipivot])
            {
                x = lhs[ipivot];
                for (pos = begin[ipivot]; pos < end[ipivot]; pos++)
                {
                    i = index[pos];
                    lhs[i] -= x * value[pos];
                    flop_count++;
                }
                if (fabs(x) > droptol)
                    pattern[nz++] = ipivot;
                else
                    lhs[ipivot] = 0.0;
            }
        }
    }
    else
    {
        for (n = 0; n < nz_symb; n++)
        {
            ipivot = pattern_symb[n];
            if (lhs[ipivot])
            {
                x = lhs[ipivot];
                for (pos = begin[ipivot]; (i = index[pos]) >= 0; pos++)
                {
                    lhs[i] -= x * value[pos];
                    flop_count++;
                }
                if (fabs(x) > droptol)
                    pattern[nz++] = ipivot;
                else
                    lhs[ipivot] = 0.0;
            }
        }
    }

    *flops += flop_count;
    return nz;
}