1 // -*- mode:C++ ; compile-command: "g++ -I.. -g -c sparse.cc" -*- 2 /* 3 * Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 4 * 5 * This program is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. 17 */ 18 #ifndef _GIAC_SPARSE_H 19 #define _GIAC_SPARSE_H 20 #include "first.h" 21 #include "index.h" 22 #include "gen.h" 23 #include <complex> 24 #include <iostream> 25 26 #ifndef NO_NAMESPACE_GIAC 27 namespace giac { 28 #endif // ndef NO_NAMESPACE_GIAC 29 // sparse matrices, not too sparse 30 struct smatrix { 31 matrice m; 32 std::vector< std::vector<int> > pos; smatrixsmatrix33 smatrix(const matrice & m_,const std::vector< std::vector<int> > & v):m(m_),pos(v) {}; smatrixsmatrix34 smatrix(){}; 35 void dbgprint() const; sizesmatrix36 int size() const { return giacmin(int(m.size()),int(pos.size())); } 37 int ncols() const; 38 }; 39 40 struct fmatrix { 41 std::vector< std::vector<giac_double> > m; 42 std::vector< std::vector<int> > pos; fmatrixfmatrix43 fmatrix(const std::vector< std::vector<giac_double> > & m_,const std::vector< std::vector<int> > & v):m(m_),pos(v) {}; fmatrixfmatrix44 fmatrix(){}; 45 void dbgprint() const; sizefmatrix46 int size() const { return giacmin(int(m.size()),int(pos.size())); } 47 int ncols() const; 48 }; 49 50 bool convert(const gen_map & d,smatrix & s); 51 bool convert(const smatrix & s,gen_map & d); 52 bool convert(const gen_map & g,vecteur & res); 53 bool convert(const vecteur & m,gen_map & res); 54 bool convert(const gen_map & d,fmatrix & s); 55 bool convert(const fmatrix & s,gen_map & d); 56 bool convert(const vecteur & source,std::vector<giac_double> & target); 57 vecteur convert(const std::vector<giac_double> & v); 58 59 void sparse_trim(const gen_map & d,gen_map &c); 60 bool need_sparse_trim(const gen_map & d); 61 62 void sparse_add(const gen_map & a,const gen_map & b,gen_map & c); 63 void sparse_neg(gen_map & c); 64 void sparse_sub(const gen_map & a,const gen_map & b,gen_map & c); 65 void sparse_mult(const gen & x,gen_map & c); 66 void sparse_div(gen_map & c,const gen & x); 67 bool sparse_mult(const gen_map & a,const gen_map & b,gen_map & c); 68 void sparse_trn(const gen_map & c,gen_map & t,bool trn,GIAC_CONTEXT); 69 void map_apply(const gen_map & a,gen_map & t,GIAC_CONTEXT,gen (* f) (const gen &,GIAC_CONTEXT) ); 70 void map_apply(const gen_map & a,const unary_function_ptr & f,gen_map & t,GIAC_CONTEXT); 71 // returns false if dimension do not match 72 bool sparse_mult(const gen_map & a,const vecteur & b,gen_map & c); 73 bool sparse_mult(const vecteur & a,const gen_map & b,gen_map & c); 74 75 bool sparse_lu(const gen_map & a,std::vector<int> & p,gen_map & l,gen_map & u_); 76 77 void sparse_mult(const smatrix & a,const vecteur & b,vecteur & c); 78 void sparse_mult(const vecteur & v,const smatrix & a,vecteur & c); 79 void sparse_mult(const std::vector<double> & v,const fmatrix & m,std::vector<double> & c); 80 void sparse_mult(const fmatrix & a,const std::vector<giac_double> & b,std::vector<giac_double> & c); 81 double l2norm(const std::vector<giac_double> & v); 82 void addvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b,std::vector<giac_double> & c); 83 void subvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b,std::vector<giac_double> & c); 84 void multvecteur(double x,const std::vector<giac_double> & a,std::vector<giac_double> & c); 85 void multvecteur(double x,std::vector<giac_double> & c); 86 std::vector<giac_double> multvecteur(double x,const std::vector<giac_double> & b); 87 std::vector<giac_double> addvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b); 88 std::vector<giac_double> subvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b); 89 90 91 // solve triangular lower inf system l*y=b 92 bool sparse_linsolve_l(const gen_map & l,const vecteur & b,vecteur & y); 93 94 // solve triangular upper system u*x=b 95 bool sparse_linsolve_u(const gen_map & u,const vecteur & b,vecteur & x); 96 97 bool is_sparse_matrix(const gen & g,int & nrows,int & ncols,int & n); 98 bool is_sparse_matrix(const gen_map & m,int & nrows,int & ncols,int & n); 99 bool is_sparse_vector(const gen & g,int & nrows,int & n); 100 bool is_sparse_vector(const gen_map & g,int & nrows,int & n); 101 102 gen sparse_conjugate_gradient(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT); 103 104 gen sparse_conjugate_gradient(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT); 105 106 gen sparse_jacobi_linsolve(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT); 107 108 gen sparse_jacobi_linsolve(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT); 109 110 std::vector<giac_double> sparse_jacobi_linsolve(const fmatrix & A,const std::vector<giac_double> & b_orig,const std::vector<giac_double> & x0,double eps,int maxiter,GIAC_CONTEXT); 111 112 std::vector<giac_double> sparse_gauss_seidel_linsolve(const fmatrix & A,const std::vector<giac_double> & b_orig,const std::vector<giac_double> & x0,double omega,double eps,int maxiter,GIAC_CONTEXT); 113 114 gen sparse_gauss_seidel_linsolve(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double omega,double eps,int maxiter,GIAC_CONTEXT); 115 116 gen sparse_gauss_seidel_linsolve(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double omega,double eps,int maxiter,GIAC_CONTEXT); 117 118 119 #ifndef NO_NAMESPACE_GIAC 120 } 121 #endif // ndef NO_NAMESPACE_GIAC 122 #endif // _GIAC_SPARSE_H 123