linear_solver.cc revision 11420:b48c0ba4f524 
1/*
2 * Copyright (c) 2015 ARM Limited
3 * All rights reserved
4 *
5 * The license below extends only to copyright in the software and shall
6 * not be construed as granting a license to any other intellectual
7 * property including but not limited to intellectual property relating
8 * to a hardware implementation of the functionality of the software
9 * licensed hereunder.  You may use the software subject to the license
10 * terms below provided that you ensure that this notice is replicated
11 * unmodified and in its entirety in all distributions of the software,
12 * modified or unmodified, in source code or in binary form.
13 *
14 * Redistribution and use in source and binary forms, with or without
15 * modification, are permitted provided that the following conditions are
16 * met: redistributions of source code must retain the above copyright
17 * notice, this list of conditions and the following disclaimer;
18 * redistributions in binary form must reproduce the above copyright
19 * notice, this list of conditions and the following disclaimer in the
20 * documentation and/or other materials provided with the distribution;
21 * neither the name of the copyright holders nor the names of its
22 * contributors may be used to endorse or promote products derived from
23 * this software without specific prior written permission.
24 *
25 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
26 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
27 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
28 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
29 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
31 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 *
37 * Authors: David Guillen Fandos
38 */
39
40#include "sim/linear_solver.hh"
41
42std::vector <double>
43LinearSystem::solve() const
44{
45    // Solve using gauss elimination, not ideal for big matrices
46    std::vector < LinearEquation > smatrix = this->matrix;
47
48    unsigned order = smatrix.size();
49    for (unsigned row = 0; row < order - 1; row++) {
50        // Look for a non-zero row, and swap
51        for (unsigned i = row; i < order; i++) {
52            if (smatrix[i][row] != 0.0f) {
53                if (i != row) {
54                    LinearEquation tmp = smatrix[i];
55                    smatrix[i] = smatrix[row];
56                    smatrix[row] = tmp;
57                }
58                break;
59            }
60        }
61
62        // Divide row by leading number to make it 1.0
63        smatrix[row] *= (1.0f / smatrix[row][row]);
64
65        // Add it (properly scaled) to the rows below
66        for (unsigned i = row + 1; i < order; i++) {
67            LinearEquation t = smatrix[row];
68            t *= -1.0f * smatrix[i][row];
69            smatrix[i] = smatrix[i] + t;
70        }
71    }
72
73    // smatrix is now a triangular matrix with diagonal being 1
74    // Just backproagate variable values from order-1 till 0
75    std::vector <double> ret(order, 0.0f);
76    for (int row = order - 1; row >= 0; row--) {
77        // Unknown value
78        ret[row] = -smatrix[row][smatrix[row].cnt()] / smatrix[row][row];
79        // Propagate variable in the cnt term
80        for (int i = row - 1; i >= 0; i--) {
81            smatrix[i][smatrix[i].cnt()] += ret[row] * smatrix[i][row];
82            smatrix[i][row] = 0.0f;
83        }
84    }
85
86    return ret;
87}
88