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
|
/****************************************************************************
* JMeshExt *
* *
* Consiglio Nazionale delle Ricerche *
* Istituto di Matematica Applicata e Tecnologie Informatiche *
* Sezione di Genova *
* IMATI-GE / CNR *
* *
* Authors: Marco Attene *
* *
* Copyright(C) 2006: IMATI-GE / CNR *
* *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
#ifndef SPARSELSYSTEM_H
#define SPARSELSYSTEM_H
#include "matrix.h"
//////////////////////////////////////////////////////////////////////////
//
// Sparse linear system
//
//////////////////////////////////////////////////////////////////////////
//! Sparse linear system Ax = B
//! Matrix A is a square matrix of size 'system_size'.
//! A is initially 0. Non-zero coefficients can be summed through
//! the method 'sumCoefficient(coeff_value, row, column)'.
//! Matrix B is a generic matrix having 'system_size' rows and 'kterm_size' columns.
//! B is initially 0. B[i][j] can be set by summing
//! values through sumKnownTerm(value, row, column).
//! A solution of the system can be obtained using one column of B trhough
//! solve(result, which_column_of_B_to_use).
//! In most geometric algorithms 'kterm_size' is 3, and the searched solution
//! is the position of some points. In this case, the value of the coordinates
//! may be retrieved by solving the system three times, once to retrieve the 'x'
//! coordinates, once for the 'y's and once for the 'z's, using the first, the
//! second and the third column of B respectively.
class sparseSystem
{
protected:
class coeffIndexPair
{
public:
int index;
double coeff;
coeffIndexPair(int a, double b) {index=a; coeff=b;}
};
class sparseSystemRow
{
public:
List cips;
sparseSystemRow() {}
~sparseSystemRow() {cips.freeNodes();}
void addCoefficient(int, double);
void print(FILE *, int);
static int rowcompare(const void *, const void *);
};
protected:
int num_equations; //!< Number of equations (rows)
int num_variables; //!< Number of variables (columns)
int kterm_size; //!< Nr. of columns of the known term
sparseSystemRow *rows; //!< Rows of the system
double **known_term; //!< Actual coefficients of the known term
public:
sparseSystem(int s, int k, int n = 0); //!< Constructs an s x n system having a k column-wide known term
~sparseSystem(); //!< Destructor
void sumCoefficient(double v, int i, int j) {rows[i].addCoefficient(j, v);} //!< Sums 'v' to A[i][j]
void setKnownTerm(double v, int i, int j) {known_term[j][i] =v;} //!< Sets j'th component of B[i] to 'v' (B[i][j]=v)
void sumKnownTerm(double v, int i, int j) {known_term[j][i]+=v;} //!< Sums j'th component of B[i] to 'v' (B[i][j]=v)
bool solve(double *solution, int j); //!< Solves the system for j'th component of B (j starts from 0). False if fails.
void print(FILE * =stdout);
};
// Sparse linear system with size-3 veriables
class sparse3System : public sparseSystem
{
bool *locks;
public:
sparse3System(int s) : sparseSystem(s, 3)
{locks = new bool[num_variables]; for (int i=0; i<num_variables; i++) locks[i]=false;}
~sparse3System() {delete locks;}
// The following calls 'solve' three times and fills the vertex array with x, y, z results
void solve(double *);
void sumKnownTerm(const double *v, int i);
void lock(int i) {locks[i]=true;}
void unlock(int i) {locks[i]=false;}
};
// Underdetermined sparse system with size-1 veriables to be solved in the least-squares sense
class leastSquaresSystem : public sparseSystem
{
bool *locks;
public:
leastSquaresSystem(int s, int n = 0);
~leastSquaresSystem() {delete locks;}
void solve(double *);
void setKnownTerm(const double v, int i) {sparseSystem::setKnownTerm(v, i, 0);}
void sumKnownTerm(const double v, int i) {sparseSystem::sumKnownTerm(v, i, 0);}
void lock(int i) {locks[i]=true;}
void unlock(int i) {locks[i]=false;}
};
#endif // SPARSELSYSTEM_H
|
