1 // -*- mode:C++ ; compile-command: "g++ -I. -I.. -I../include -g -c sparse.cc -fno-strict-aliasing -DGIAC_GENERIC_CONSTANTS -DHAVE_CONFIG_H -DIN_GIAC" -*- 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2000,14 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 using namespace std; 20 #include <cmath> 21 #include <stdexcept> 22 #include <map> 23 #include <iostream> 24 #include "gen.h" 25 #include "sparse.h" 26 #include "vecteur.h" 27 #include "modpoly.h" 28 #include "unary.h" 29 #include "symbolic.h" 30 #include "usual.h" 31 #include "sym2poly.h" 32 #include "solve.h" 33 #include "prog.h" 34 #include "subst.h" 35 #include "permu.h" 36 #include "plot.h" 37 #include "misc.h" 38 #include "ti89.h" 39 #include "csturm.h" 40 #include "giacintl.h" 41 #ifdef HAVE_LIBGSL 42 #include <gsl/gsl_linalg.h> 43 #include <gsl/gsl_eigen.h> 44 #include <gsl/gsl_poly.h> 45 #endif 46 47 #ifndef NO_NAMESPACE_GIAC 48 namespace giac { 49 #endif // ndef NO_NAMESPACE_GIAC 50 sparse_add(const gen_map & a,const gen_map & b,gen_map & c)51 void sparse_add(const gen_map & a,const gen_map & b,gen_map & c){ 52 c.clear(); 53 comparegen key; 54 gen_map::const_iterator it=a.begin(),itend=a.end(),jt=b.begin(),jtend=b.end(); 55 for (;it!=itend && jt!=jtend;){ 56 if (it->first==jt->first){ 57 gen tmp=it->second+jt->second; 58 if (!is_zero(tmp)){ 59 c[it->first]=tmp; 60 } 61 ++it; 62 ++jt; 63 continue; 64 } 65 if (key(it->first,jt->first)){ 66 c[it->first]=it->second; 67 ++it; 68 continue; 69 } 70 c[jt->first]=jt->second; 71 ++jt; 72 } 73 for (;it!=itend;++it){ 74 c[it->first]=it->second; 75 } 76 for (;jt!=jtend;++jt){ 77 c[jt->first]=jt->second; 78 } 79 } 80 sparse_sub(const gen_map & a,const gen_map & b,gen_map & c)81 void sparse_sub(const gen_map & a,const gen_map & b,gen_map & c){ 82 c.clear(); 83 comparegen key; 84 gen_map::const_iterator it=a.begin(),itend=a.end(),jt=b.begin(),jtend=b.end(); 85 for (;it!=itend && jt!=jtend;){ 86 if (it->first==jt->first){ 87 gen tmp=it->second-jt->second; 88 if (!is_zero(tmp)){ 89 c[it->first]=tmp; 90 } 91 ++it; 92 ++jt; 93 continue; 94 } 95 if (key(it->first,jt->first)){ 96 c[it->first]=it->second; 97 ++it; 98 continue; 99 } 100 c[jt->first]=-jt->second; 101 ++jt; 102 } 103 for (;it!=itend;++it){ 104 c[it->first]=it->second; 105 } 106 for (;jt!=jtend;++jt){ 107 c[jt->first]=-jt->second; 108 } 109 } 110 sparse_mult(const gen & x,gen_map & c)111 void sparse_mult(const gen & x,gen_map & c){ 112 if (is_zero(x)){ 113 c.clear(); 114 return; 115 } 116 gen_map::iterator it=c.begin(),itend=c.end(); 117 for (;it!=itend;++it){ 118 it->second = x*it->second; 119 } 120 } 121 sparse_div(gen_map & c,const gen & x)122 void sparse_div(gen_map & c,const gen & x){ 123 if (is_inf(x)){ 124 c.clear(); 125 return; 126 } 127 gen_map::iterator it=c.begin(),itend=c.end(); 128 for (;it!=itend;++it){ 129 it->second = it->second/x; 130 } 131 } 132 sparse_neg(gen_map & c)133 void sparse_neg(gen_map & c){ 134 gen_map::iterator it=c.begin(),itend=c.end(); 135 for (;it!=itend;++it){ 136 it->second = -it->second; 137 } 138 } 139 140 struct sparse_line_begin_t { 141 longlong i; 142 gen_map::const_iterator begin,end; sparse_line_begin_tgiac::sparse_line_begin_t143 sparse_line_begin_t(longlong i_,gen_map::const_iterator it_,gen_map::const_iterator jt_):i(i_),begin(it_),end(jt_){}; sparse_line_begin_tgiac::sparse_line_begin_t144 sparse_line_begin_t(){}; 145 }; 146 operator <(const sparse_line_begin_t & a,const sparse_line_begin_t & b)147 bool operator < (const sparse_line_begin_t & a,const sparse_line_begin_t & b){ 148 return a.i<b.i; 149 } 150 dicho(const std::vector<sparse_line_begin_t> & B,longlong pos,gen_map::const_iterator & begin,gen_map::const_iterator & end)151 bool dicho(const std::vector<sparse_line_begin_t> & B,longlong pos,gen_map::const_iterator & begin,gen_map::const_iterator & end){ 152 if (B.empty()) 153 return false; 154 longlong a=0,b=B.size()-1,c; // we are between a and b 155 if (pos<B[0].i || pos>B[b].i) 156 return false; 157 for (;a<b-1;){ 158 c=(a+b)/2; 159 if (pos<B[c].i) 160 b=c; 161 else 162 a=c; 163 } 164 // a>=b-1 hence a==b-1 or a==b 165 if (pos==B[a].i){ 166 begin=B[a].begin; 167 end=B[a].end; 168 for (;begin!=end;++begin){ 169 if (begin->second!=0) 170 break; 171 } 172 return true; 173 } 174 if (pos==B[b].i){ 175 begin=B[b].begin; 176 end=B[b].end; 177 for (;begin!=end;++begin){ 178 if (begin->second!=0) 179 break; 180 } 181 return true; 182 } 183 return false; 184 } 185 sparse_trim(const gen_map & d,gen_map & c)186 void sparse_trim(const gen_map & d,gen_map &c){ 187 gen_map::const_iterator it=d.begin(),itend=d.end(); 188 for (;it!=itend;++it){ 189 if (!is_zero(it->second)) 190 c[it->first]=it->second; 191 } 192 } 193 need_sparse_trim(const gen_map & d)194 bool need_sparse_trim(const gen_map & d){ 195 gen_map::const_iterator it=d.begin(),itend=d.end(); 196 for (;it!=itend;++it){ 197 if (is_zero(it->second)) 198 return true; 199 } 200 return false; 201 } 202 find_line_begin(gen_map::const_iterator jt,gen_map::const_iterator jtend,std::vector<sparse_line_begin_t> & B)203 void find_line_begin(gen_map::const_iterator jt,gen_map::const_iterator jtend,std::vector<sparse_line_begin_t> & B){ 204 B.clear(); 205 gen_map::const_iterator begin; 206 longlong prev=-1; 207 for (;jt!=jtend;++jt){ 208 longlong cur=jt->first._VECTptr->front().val; 209 if (cur==prev) 210 continue; 211 if (prev!=-1 && begin!=jtend) 212 B.push_back(sparse_line_begin_t(prev,begin,jt)); 213 prev=cur; 214 begin=jt; 215 } 216 if (begin!=jtend) 217 B.push_back(sparse_line_begin_t(prev,begin,jt)); 218 } 219 sparse_mult(const gen_map & a,const gen_map & b,gen_map & c)220 bool sparse_mult(const gen_map & a,const gen_map & b,gen_map & c){ 221 c.clear(); 222 if (a.empty() || b.empty()) 223 return true; 224 gen_map::const_iterator it=a.begin(),itend=a.end(),jt=b.begin(),jtend=b.end(),begin,end; 225 if (jt->first.type==_INT_){ 226 gen_map d; 227 // matrix * std::vector 228 for (;it!=itend;++it){ 229 jt=b.find(it->first._VECTptr->back()); 230 if (jt!=jtend) 231 d[it->first._VECTptr->front()] += it->second*jt->second; 232 } 233 sparse_trim(d,c); 234 return true; 235 } 236 if (it->first.type==_INT_){ 237 gen_map d; 238 std::vector<sparse_line_begin_t> B; 239 // compute positions of line begin of b in B 240 find_line_begin(jt,jtend,B); 241 for (;it!=itend;++it){ 242 if (!dicho(B,it->first.val,begin,end)) 243 continue; 244 gen coeffa=it->second; 245 for (;begin!=end;++begin){ 246 d[begin->first._VECTptr->back()] += coeffa*begin->second; 247 } 248 } 249 sparse_trim(d,c); 250 return true; 251 } 252 if (1){ 253 smatrix B; 254 if (!convert(b,B)) 255 return false; 256 int C=B.ncols(); 257 int previ=-1; 258 vecteur cur(C); 259 for (;;){ 260 int i=-2,k=-1; 261 gen coeffa; 262 if (it!=itend){ 263 i=it->first._VECTptr->front().val; 264 k=it->first._VECTptr->back().val; 265 coeffa=it->second; 266 ++it; 267 } 268 if (i!=previ){ 269 if (previ>=0){ 270 // copy cur in c 271 vecteur pos(2,previ); 272 for (int j=0;j<C;++j){ 273 gen & cj=cur[j]; 274 if (is_zero(cj)) 275 continue; 276 pos[1]=j; 277 c[gen(pos,_SEQ__VECT)]=cj; 278 cj=0; 279 } 280 } 281 if (i==-2 && it==itend) break; 282 previ=i; 283 } 284 if (k>=B.m.size()) 285 continue; 286 // a_[i,k]*B[k,j] 287 const_iterateur bt=B.m[k]._VECTptr->begin(),btend=B.m[k]._VECTptr->end(); 288 std::vector<int>::const_iterator bpos=B.pos[k].begin(); 289 for (;bt!=btend;++bpos,++bt){ 290 cur[*bpos] += coeffa*(*bt); 291 } 292 } // end for 293 } 294 else { 295 gen_map d; 296 std::vector<sparse_line_begin_t> B; 297 // compute positions of line begin of b in B 298 find_line_begin(jt,jtend,B); 299 // sort(B.begin(),B.end()); 300 vecteur pos(2); 301 for (;it!=itend;++it){ 302 pos[0]=it->first._VECTptr->front(); 303 longlong cur=it->first._VECTptr->back().val; 304 gen coeffa=it->second; 305 // dichotomic search in B 306 if (!dicho(B,cur,begin,end)) 307 continue; 308 for (;begin!=end;++begin){ 309 pos[1]=begin->first._VECTptr->back(); 310 d[gen(pos,_SEQ__VECT)] += coeffa*begin->second; 311 } 312 } 313 sparse_trim(d,c); 314 } 315 return true; 316 } 317 sparse_mult(const gen_map & a,const vecteur & b,gen_map & c)318 bool sparse_mult(const gen_map & a,const vecteur & b,gen_map & c){ 319 c.clear(); 320 if (a.empty() || b.empty()) 321 return true; 322 gen_map d; 323 gen_map::const_iterator it=a.begin(),itend=a.end(); 324 vecteur pos(2); 325 for (;it!=itend;++it){ 326 gen i=it->first._VECTptr->front(); 327 int j=it->first._VECTptr->back().val; 328 gen coeffa=it->second; 329 if (j>=b.size()) 330 return false; 331 gen bj=b[j]; 332 if (bj.type!=_VECT) 333 d[i] += coeffa*bj; 334 else { 335 pos[0]=i; 336 const_iterateur jt=bj._VECTptr->begin(),jtend=bj._VECTptr->end(); 337 for (int k=0;jt!=jtend;++jt,++k){ 338 pos[1]=k; 339 d[gen(pos,_SEQ__VECT)] += coeffa*(*jt); 340 } 341 } 342 } 343 sparse_trim(d,c); 344 return true; 345 } 346 sparse_mult(const smatrix & a,const vecteur & b,vecteur & c)347 void sparse_mult(const smatrix & a,const vecteur & b,vecteur & c){ 348 c.clear(); 349 int l=a.size(); 350 c.reserve(l); 351 for (int i=0;i<l;++i){ 352 gen tmp=0; 353 const_iterateur at=a.m[i]._VECTptr->begin(),atend=a.m[i]._VECTptr->end(); 354 std::vector<int>::const_iterator apos=a.pos[i].begin(); 355 for (;at!=atend;++apos,++at){ 356 tmp += (*at)*b[*apos]; 357 } 358 c.push_back(tmp); 359 } 360 } 361 sparse_mult(const fmatrix & a,const std::vector<giac_double> & b,std::vector<giac_double> & c)362 void sparse_mult(const fmatrix & a,const std::vector<giac_double> & b,std::vector<giac_double> & c){ 363 c.clear(); 364 int l=a.size(); 365 c.reserve(l); 366 for (int i=0;i<l;++i){ 367 double tmp=0; 368 std::vector<double>::const_iterator at=a.m[i].begin(),atend=a.m[i].end(); 369 std::vector<int>::const_iterator apos=a.pos[i].begin(); 370 for (;at!=atend;++apos,++at){ 371 tmp += (*at)*b[*apos]; 372 } 373 c.push_back(tmp); 374 } 375 } 376 sparse_mult(const vecteur & v,const smatrix & m,vecteur & c)377 void sparse_mult(const vecteur & v,const smatrix & m,vecteur & c){ 378 c.clear(); 379 int l=m.size(); 380 c.resize(l); 381 for (int i=0;i<l;++i){ 382 gen coeff=v[i]; 383 const_iterateur mt=m.m[i]._VECTptr->begin(),mtend=m.m[i]._VECTptr->end(); 384 std::vector<int>::const_iterator mpos=m.pos[i].begin(); 385 for (;mt!=mtend;++mpos,++mt){ 386 c[*mpos] += coeff*(*mt); 387 } 388 } 389 } 390 sparse_mult(const std::vector<double> & v,const fmatrix & m,std::vector<double> & c)391 void sparse_mult(const std::vector<double> & v,const fmatrix & m,std::vector<double> & c){ 392 c.clear(); 393 int l=m.size(); 394 c.resize(l); 395 for (int i=0;i<l;++i){ 396 double coeff=v[i]; 397 std::vector<double>::const_iterator mt=m.m[i].begin(),mtend=m.m[i].end(); 398 std::vector<int>::const_iterator mpos=m.pos[i].begin(); 399 for (;mt!=mtend;++mpos,++mt){ 400 c[*mpos] += coeff*(*mt); 401 } 402 } 403 } 404 sparse_mult(const vecteur & a,const gen_map & b,gen_map & c)405 bool sparse_mult(const vecteur & a,const gen_map & b,gen_map & c){ 406 c.clear(); 407 if (a.empty() || b.empty()) 408 return true; 409 bool mat=ckmatrix(a); vecteur a_; 410 if (mat) 411 a_=mtran(a); 412 gen_map d; 413 gen_map::const_iterator it=b.begin(),itend=b.end(); 414 vecteur pos(2); 415 for (;it!=itend;++it){ 416 int j=it->first._VECTptr->front().val; 417 gen k=it->first._VECTptr->back(); 418 gen coeffa=it->second; 419 if (!mat){ 420 if (j>=a.size()) 421 return false; 422 d[k] += a[j]*coeffa; 423 } 424 else { 425 if (j>=a_.size()) 426 return false; 427 pos[1]=k; 428 const_iterateur jt=a_[j]._VECTptr->begin(),jtend=a_[j]._VECTptr->end(); 429 for (int i=0;jt!=jtend;++jt,++i){ 430 pos[0]=i; 431 d[gen(pos,_SEQ__VECT)] += (*jt)*coeffa; 432 } 433 } 434 } 435 sparse_trim(d,c); 436 return true; 437 } 438 439 // transpose or transconjugate sparse_trn(const gen_map & c,gen_map & t,bool trn,GIAC_CONTEXT)440 void sparse_trn(const gen_map & c,gen_map & t,bool trn,GIAC_CONTEXT){ 441 t.clear(); 442 gen_map::const_iterator it=c.begin(),itend=c.end(); 443 for (;it!=itend;++it){ 444 gen i=it->first; 445 if (i.type==_INT_) 446 i=makesequence(0,i); 447 else 448 i=makesequence(i._VECTptr->back(),i._VECTptr->front()); 449 t[i]= trn?conj(it->second,contextptr):it->second; 450 } 451 } 452 map_apply(const gen_map & c,gen_map & t,GIAC_CONTEXT,gen (* f)(const gen &,GIAC_CONTEXT))453 void map_apply(const gen_map & c,gen_map & t,GIAC_CONTEXT,gen (* f) (const gen &,GIAC_CONTEXT) ){ 454 t.clear(); 455 gen_map::const_iterator it=c.begin(),itend=c.end(); 456 for (;it!=itend;++it){ 457 gen g=f(it->second,contextptr); 458 if (!is_zero(g)) 459 t[it->first]=g; 460 } 461 } 462 map_apply(const gen_map & c,const unary_function_ptr & f,gen_map & t,GIAC_CONTEXT)463 void map_apply(const gen_map & c,const unary_function_ptr & f,gen_map & t,GIAC_CONTEXT){ 464 t.clear(); 465 gen_map::const_iterator it=c.begin(),itend=c.end(); 466 for (;it!=itend;++it){ 467 gen g=f(it->second,contextptr); 468 if (!is_zero(g)) 469 t[it->first]=g; 470 } 471 } 472 sparse_swaprows(gen_map & u,std::vector<sparse_line_begin_t> & B,int r1,int r2,int c=-1)473 bool sparse_swaprows(gen_map & u,std::vector<sparse_line_begin_t>& B,int r1,int r2,int c=-1){ 474 if (r1>r2) 475 return sparse_swaprows(u,B,r2,r1,c); 476 if (r1==r2) 477 return true; 478 gen_map::const_iterator r1b0,r1b,r1e,r2b0,r2b,r2e; 479 bool dicho1=dicho(B,r1,r1b0,r1e),dicho2=dicho(B,r2,r2b0,r2e); 480 if (!dicho1 && !dicho2) 481 return false; 482 vecteur c1; 483 if (dicho1){ 484 for (r1b=r1b0;r1b!=r1e;++r1b){ 485 gen & C=r1b->first._VECTptr->back(); 486 if (c>=0 && C.val>c){ 487 r1e=r1b; 488 break; 489 } 490 c1.push_back(C); 491 c1.push_back(r1b->second); 492 } 493 } 494 vecteur c2; 495 if (dicho2){ 496 for (r2b=r2b0;r2b!=r2e;++r2b){ 497 gen & C=r2b->first._VECTptr->back(); 498 if (c>=0 && C.val>c){ 499 r2e=r2b; 500 break; 501 } 502 c2.push_back(C); 503 c2.push_back(r2b->second); 504 } 505 } 506 if (dicho1) 507 u.erase(*(gen_map::iterator *) &r1b0, *(gen_map::iterator *) &r1e); 508 if (dicho2) 509 u.erase(*(gen_map::iterator *) &r2b0, *(gen_map::iterator *) &r2e); 510 for (int i=0;i<int(c2.size());i+=2){ 511 u[makesequence(r1,c2[i])]=c2[i+1]; 512 } 513 for (int i=0;i<int(c1.size());i+=2){ 514 u[makesequence(r2,c1[i])]=c1[i+1]; 515 } 516 return true; 517 } 518 sparse_rowadd(gen_map & u,int r,int col,gen_map::const_iterator it,gen_map::const_iterator itend,const gen & coeff,gen_map::const_iterator jt,gen_map::const_iterator jtend)519 void sparse_rowadd(gen_map & u,int r,int col,gen_map::const_iterator it,gen_map::const_iterator itend,const gen & coeff,gen_map::const_iterator jt,gen_map::const_iterator jtend){ 520 vecteur c; 521 gen_map::const_iterator itsave=it; 522 for (;it!=itend && jt!=jtend;){ 523 int ic=it->first._VECTptr->back().val,jc=jt->first._VECTptr->back().val; 524 if (ic<col){ 525 ++it; 526 continue; 527 } 528 if (jc<col){ 529 ++it; 530 continue; 531 } 532 if (ic==jc){ 533 gen tmp=it->second+coeff*jt->second; 534 c.push_back(ic); 535 c.push_back(tmp); 536 ++it; 537 ++jt; 538 continue; 539 } 540 if (ic<jc){ 541 c.push_back(ic); 542 c.push_back(it->second); 543 ++it; 544 continue; 545 } 546 c.push_back(jc); 547 c.push_back(coeff*jt->second); 548 ++jt; 549 } 550 for (;it!=itend;++it){ 551 int ic=it->first._VECTptr->back().val; 552 if (ic<col) 553 continue; 554 c.push_back(ic); 555 c.push_back(it->second); 556 } 557 for (;jt!=jtend;++jt){ 558 int jc=jt->first._VECTptr->back().val; 559 if (jc<col) 560 continue; 561 c.push_back(jc); 562 c.push_back(coeff*jt->second); 563 } 564 // copy c in the map 565 // should erase 0 566 for (unsigned i=0;i<c.size();i+=2){ 567 u[makesequence(r,c[i])]=c[i+1]; 568 } 569 } 570 sparse_lu(const gen_map & a,std::vector<int> & p,gen_map & l,gen_map & u_)571 bool sparse_lu(const gen_map & a,std::vector<int> & p,gen_map & l,gen_map & u_){ 572 int L,C,n; 573 if (!is_sparse_matrix(a,L,C,n)) 574 return false; 575 l.clear(); 576 gen_map u(a); 577 p=std::vector<int>(L); 578 for (int i=0;i<L;++i) 579 p[i]=i; 580 int r=0,c=0; 581 for (;r<L && c<C;){ 582 std::vector<sparse_line_begin_t> B; 583 find_line_begin(u.begin(),u.end(),B); 584 // search pivot in column c starting at line l 585 gen_map::const_iterator begin,end; 586 int lp=r,pivotline=-1; 587 gen pivot=0; 588 for (;lp<L;++lp){ 589 if (!dicho(B,lp,begin,end)) 590 continue; 591 for (;begin!=end;++begin){ 592 if (begin->first._VECTptr->back()==c){ 593 gen temp=begin->second; 594 if (temp.islesscomplexthan(pivot)){ 595 pivot=temp; 596 pivotline=lp; 597 } 598 break; 599 } 600 if (begin->first._VECTptr->back()>c) 601 break; 602 } 603 } 604 if (pivotline==-1){ 605 ++c; 606 continue; 607 } 608 // exchange rows 609 if (pivotline!=r){ 610 //CERR << "pivotline,r " << pivotline << "," << r << " col " << c << '\n'; 611 //CERR << "avant " << l << '\n' << u << '\n'; 612 sparse_swaprows(u,B,r,pivotline); 613 find_line_begin(l.begin(),l.end(),B); 614 sparse_swaprows(l,B,r,pivotline,c); 615 find_line_begin(u.begin(),u.end(),B); 616 //CERR << "apres " << l << '\n' << u << '\n'; 617 swapint(p[r],p[pivotline]); 618 } 619 // reduce rows and fill l 620 l[makesequence(r,c)]=1; 621 gen_map::const_iterator begin2,end2; 622 if (!dicho(B,r,begin2,end2)) 623 return false; 624 for (int R=r+1;R<L;++R){ 625 if (!dicho(B,R,begin,end)) 626 continue; 627 if (begin->first._VECTptr->back()!=c) 628 continue; 629 gen coeff=begin->second/pivot; 630 l[makesequence(R,c)]=coeff; 631 sparse_rowadd(u,R,c,begin,end,-coeff,begin2,end2); 632 } 633 ++r; ++c; 634 } // end r<L && c<C 635 sparse_trim(u,u_); 636 return true; 637 } 638 639 // solve triangular lower inf system l*y=b sparse_linsolve_l(const gen_map & l,const vecteur & b,vecteur & y)640 bool sparse_linsolve_l(const gen_map & l,const vecteur & b,vecteur & y){ 641 int n=int(b.size()); 642 y.resize(n); 643 std::vector<sparse_line_begin_t> L; 644 find_line_begin(l.begin(),l.end(),L); 645 gen_map::const_iterator begin,end; 646 for (int i=0;i<n;++i){ 647 if (!dicho(L,i,begin,end)) 648 return false; 649 gen res=b[i]; bool ok=false; 650 for (;begin!=end;++begin){ 651 int j=begin->first._VECTptr->back().val; 652 if (j>i) 653 return false; 654 if (j==i){ 655 res=res/begin->second; 656 ok=true; 657 } 658 else 659 res -= y[j]*begin->second; 660 } 661 if (!ok) 662 return false; 663 y[i]=res; 664 } 665 return true; 666 } 667 sparse_linsolve_l(const fmatrix & l,const std::vector<giac_double> & b,std::vector<giac_double> & y)668 bool sparse_linsolve_l(const fmatrix & l,const std::vector<giac_double> & b,std::vector<giac_double> & y){ 669 int n=int(b.size()); 670 y.resize(n); 671 for (int i=0;i<n;++i){ 672 const std::vector<int> & posi=l.pos[i]; 673 const std::vector<giac_double> & li = l.m[i]; 674 double res=b[i]; int s=int(posi.size()); bool ok=false; 675 for (int j=0;j<s;++j){ 676 int pos=posi[j]; 677 if (pos>i) 678 return false; 679 if (pos==i){ 680 res /= li[j]; 681 ok=true; 682 } 683 else 684 res -= y[pos]*li[j]; 685 } 686 if (!ok) 687 return false; 688 y[i]=res; 689 } 690 return true; 691 } 692 693 // solve triangular upper system u*x=b sparse_linsolve_u(const gen_map & u,const vecteur & b,vecteur & x)694 bool sparse_linsolve_u(const gen_map & u,const vecteur & b,vecteur & x){ 695 int n=int(b.size()); 696 x.resize(n); 697 std::vector<sparse_line_begin_t> U; 698 find_line_begin(u.begin(),u.end(),U); 699 gen_map::const_iterator begin,end; 700 for (int i=n-1;i>=0;--i){ 701 if (!dicho(U,i,begin,end)) 702 return false; 703 gen res=b[i],coeff; bool ok=false; 704 for (;begin!=end;++begin){ 705 int j=begin->first._VECTptr->back().val; 706 if (j<i) 707 return false; 708 if (j==i){ 709 coeff=begin->second; 710 ok=true; 711 } 712 else 713 res -= x[j]*begin->second; 714 } 715 if (!ok) 716 return false; 717 x[i]=res/coeff; 718 } 719 return true; 720 } 721 dbgprint() const722 void smatrix::dbgprint() const { 723 gen_map d; convert(*this,d); 724 CERR << d << '\n'; 725 } 726 dbgprint() const727 void fmatrix::dbgprint() const { 728 for (int i=0;i<int(pos.size());++i){ 729 const std::vector<int> & posi=pos[i]; 730 const std::vector<double> & mi=m[i]; 731 CERR << "line " << i << ": "; 732 for (int j=0;j<posi.size();++j){ 733 CERR << posi[j]<<"="<<mi[j] <<", "; 734 } 735 CERR << '\n'; 736 } 737 } 738 convert(const gen_map & d,smatrix & s)739 bool convert(const gen_map & d,smatrix & s){ 740 int nrows,ncols,n; 741 if (!is_sparse_matrix(d,nrows,ncols,n)) 742 return false; 743 s.pos=std::vector< std::vector<int> >(nrows); 744 s.m=vecteur(nrows); 745 for (int i=0;i<nrows;++i) 746 s.m[i]=vecteur(0); 747 int prev=-1; 748 gen_map::const_iterator it=d.begin(),itend=d.end(); 749 for (;it!=itend;++it){ 750 gen p=it->first; 751 if (p.type!=_VECT || p._VECTptr->size()!=2) 752 return false; 753 int l=p._VECTptr->front().val,c=p._VECTptr->back().val; 754 s.pos[l].push_back(c); 755 s.m[l]._VECTptr->push_back(it->second); 756 } 757 return true; 758 } 759 convert(const gen_map & d,fmatrix & s)760 bool convert(const gen_map & d,fmatrix & s){ 761 int nrows,ncols,n; 762 if (!is_sparse_matrix(d,nrows,ncols,n)) 763 return false; 764 s.pos=std::vector< std::vector<int> >(nrows); 765 s.m=std::vector< std::vector<giac_double> >(nrows); 766 int prev=-1; 767 gen_map::const_iterator it=d.begin(),itend=d.end(); 768 for (;it!=itend;++it){ 769 gen p=it->first; 770 if (p.type!=_VECT || p._VECTptr->size()!=2) 771 return false; 772 int l=p._VECTptr->front().val,c=p._VECTptr->back().val; 773 s.pos[l].push_back(c); 774 gen tmp=evalf_double(it->second,1,context0); 775 if (tmp.type!=_DOUBLE_) 776 return false; 777 s.m[l].push_back(tmp._DOUBLE_val); 778 } 779 return true; 780 } 781 convert(const smatrix & s,gen_map & d)782 bool convert(const smatrix & s,gen_map & d){ 783 d.clear(); 784 int ss=s.size(); 785 for (int i=0;i<ss;++i){ 786 const vecteur & v =*s.m[i]._VECTptr; 787 const std::vector<int> & p=s.pos[i]; 788 if (v.size()!=p.size()) 789 return false; 790 const_iterateur it=v.begin(),itend=v.end(); 791 std::vector<int>::const_iterator jt=p.begin(); 792 for (;it!=itend;++jt,++it){ 793 d[makesequence(i,*jt)]=*it; 794 } 795 } 796 return true; 797 } 798 convert(const fmatrix & s,gen_map & d)799 bool convert(const fmatrix & s,gen_map & d){ 800 d.clear(); 801 int ss=s.size(); 802 for (int i=0;i<ss;++i){ 803 const std::vector<double> & v =s.m[i]; 804 const std::vector<int> & p=s.pos[i]; 805 if (v.size()!=p.size()) 806 return false; 807 std::vector<double>::const_iterator it=v.begin(),itend=v.end(); 808 std::vector<int>::const_iterator jt=p.begin(); 809 for (;it!=itend;++jt,++it){ 810 d[makesequence(i,*jt)]=*it; 811 } 812 } 813 return true; 814 } 815 convert(const std::vector<giac_double> & v)816 vecteur convert(const std::vector<giac_double> & v){ 817 vecteur res; 818 res.reserve(v.size()); 819 for (int i=0;i<v.size();++i) 820 res.push_back(v[i]); 821 return res; 822 } 823 convert(const vecteur & source,std::vector<giac_double> & target)824 bool convert(const vecteur & source,std::vector<giac_double> & target){ 825 const_iterateur it=source.begin(),itend=source.end(); 826 target.clear(); 827 target.reserve(itend-it); 828 for (;it!=itend;++it){ 829 gen tmp=evalf_double(*it,1,context0); 830 if (tmp.type!=_DOUBLE_) 831 return false; 832 target.push_back(tmp._DOUBLE_val); 833 } 834 return true; 835 } 836 convert(const vecteur & m,gen_map & res)837 bool convert(const vecteur & m,gen_map & res){ 838 if (ckmatrix(m)){ 839 for (int i=0;i<int(m.size());++i){ 840 const vecteur & v = *m[i]._VECTptr; 841 for (int j=0;j<int(v.size());++j){ 842 if (!is_zero(v[j])) 843 res[makesequence(i,j)]=v[j]; 844 } 845 } 846 } 847 else { 848 for (int i=0;i<int(m.size());++i){ 849 if (!is_zero(m[i])) 850 res[i]=m[i]; 851 } 852 } 853 return true; 854 } 855 convert(const gen_map & g,vecteur & res)856 bool convert(const gen_map & g,vecteur & res){ 857 int n,nrows,ncols; 858 if (!is_sparse_matrix(g,nrows,ncols,n)){ 859 if (!is_sparse_vector(g,nrows,n)) 860 return false; 861 res=vecteur(nrows); 862 gen_map::const_iterator it=g.begin(),itend=g.end(); 863 for (;it!=itend;++it){ 864 gen l=it->first; 865 is_integral(l); 866 res[l.val]=it->second; 867 } 868 return true; 869 } 870 res=vecteur(nrows); 871 for (int i=0;i<nrows;++i){ 872 vecteur l(ncols); 873 res[i]=l; 874 } 875 gen_map::const_iterator it=g.begin(),itend=g.end(); 876 for (;it!=itend;++it){ 877 gen G=it->first; 878 gen l=G._VECTptr->front(); 879 gen c=G._VECTptr->back(); 880 is_integral(l); is_integral(c); 881 (*res[l.val]._VECTptr)[c.val]=it->second; 882 } 883 return true; 884 } 885 is_sparse_matrix(const gen & g,int & nrows,int & ncols,int & n)886 bool is_sparse_matrix(const gen & g,int & nrows,int & ncols,int & n){ 887 if (g.type!=_MAP) 888 return false; 889 gen_map & m=*g._MAPptr; 890 return is_sparse_matrix(*g._MAPptr,nrows,ncols,n); 891 } is_sparse_matrix(const gen_map & m,int & nrows,int & ncols,int & n)892 bool is_sparse_matrix(const gen_map & m,int & nrows,int & ncols,int & n){ 893 nrows=0;ncols=0;n=0; 894 gen_map::const_iterator it=m.begin(),itend=m.end(); 895 for (;it!=itend;++n,++it){ 896 gen g=it->first; 897 if (g.type!=_VECT || g._VECTptr->size()!=2) 898 return false; 899 gen l=g._VECTptr->front(); 900 gen c=g._VECTptr->back(); 901 if (!is_integral(l) || !is_integral(c) || l.val<0 || c.val<0) 902 return false; 903 if (nrows<=l.val) 904 nrows=l.val+1; 905 if (ncols<=c.val) 906 ncols=c.val+1; 907 } 908 return true; 909 } 910 is_sparse_vector(const gen & g,int & nrows,int & n)911 bool is_sparse_vector(const gen & g,int & nrows,int & n){ 912 if (g.type!=_MAP) 913 return false; 914 return is_sparse_vector(*g._MAPptr,nrows,n); 915 } is_sparse_vector(const gen_map & m,int & nrows,int & n)916 bool is_sparse_vector(const gen_map & m,int & nrows,int & n){ 917 nrows=0;n=0; 918 gen_map::const_iterator it=m.begin(),itend=m.end(); 919 for (;it!=itend;++n,++it){ 920 gen l=it->first; 921 if (!is_integral(l) || l.val<0) 922 return false; 923 if (nrows<=l.val) 924 nrows=l.val+1; 925 } 926 return true; 927 } 928 ncols_(const std::vector<std::vector<int>> & pos)929 static int ncols_(const std::vector< std::vector<int> > & pos){ 930 int res=-1; 931 for (unsigned i=0;i<pos.size();++i){ 932 const std::vector<int> & pi=pos[i]; 933 if (!pi.empty()) 934 res=giacmax(res,pi.back()); 935 } 936 return res+1; 937 } 938 ncols() const939 int smatrix::ncols() const { 940 return ncols_(pos); 941 } 942 ncols() const943 int fmatrix::ncols() const { 944 return ncols_(pos); 945 } 946 sparse_conjugate_gradient(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT)947 gen sparse_conjugate_gradient(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT){ 948 int n=int(b_orig.size()); 949 vecteur tmp(n); 950 sparse_mult(A,x0,tmp); 951 vecteur b=subvecteur(b_orig,tmp); 952 vecteur xk(x0); 953 vecteur rk(b),pk(b); 954 gen rk2=scalarproduct(rk,rk,contextptr); 955 vecteur Apk(n); 956 for (int k=1;k<=maxiter;++k){ 957 sparse_mult(A,pk,Apk); 958 gen alphak=rk2/scalarproduct(pk,Apk,contextptr); 959 multvecteur(alphak,pk,tmp); 960 addvecteur(xk,tmp,xk); 961 multvecteur(alphak,Apk,tmp); 962 subvecteur(rk,tmp,rk); 963 gen newrk2=scalarproduct(rk,rk,contextptr); 964 if (is_greater(eps*eps,newrk2,contextptr)) 965 return xk; 966 multvecteur(newrk2/rk2,pk,tmp); 967 addvecteur(rk,tmp,pk); 968 rk2=newrk2; 969 } 970 *logptr(contextptr) << gettext("Warning! Leaving conjugate gradient algorithm after dimension of matrix iterations. Check that your matrix is hermitian/symmetric definite.") << '\n'; 971 return xk; 972 } 973 sparse_conjugate_gradient(const fmatrix & A,const std::vector<giac_double> & b_orig,const std::vector<giac_double> & x0,double eps,int maxiter,GIAC_CONTEXT)974 std::vector<giac_double> sparse_conjugate_gradient(const fmatrix & A,const std::vector<giac_double> & b_orig,const std::vector<giac_double> & x0,double eps,int maxiter,GIAC_CONTEXT){ 975 int n=int(b_orig.size()); 976 std::vector<giac_double> tmp(n); 977 sparse_mult(A,x0,tmp); 978 std::vector<giac_double> b; 979 subvecteur(b_orig,tmp,b); 980 std::vector<giac_double> xk(x0); 981 std::vector<giac_double> rk(b),pk(b); 982 giac_double rk2=dotvecteur(rk,rk); 983 std::vector<giac_double> Apk(n); 984 for (int k=1;k<=maxiter;++k){ 985 sparse_mult(A,pk,Apk); 986 giac_double alphak=rk2/dotvecteur(pk,Apk); 987 multvecteur(alphak,pk,tmp); 988 addvecteur(xk,tmp,xk); 989 multvecteur(alphak,Apk,tmp); 990 subvecteur(rk,tmp,rk); 991 giac_double newrk2=dotvecteur(rk,rk); 992 if (eps*eps>=newrk2) 993 return xk; 994 multvecteur(newrk2/rk2,pk,tmp); 995 addvecteur(rk,tmp,pk); 996 rk2=newrk2; 997 } 998 *logptr(contextptr) << gettext("Warning! Leaving conjugate gradient algorithm after dimension of matrix iterations. Check that your matrix is hermitian/symmetric definite.") << '\n'; 999 return xk; 1000 } 1001 sparse_conjugate_gradient(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT)1002 gen sparse_conjugate_gradient(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT){ 1003 if (has_num_coeff(b_orig)){ 1004 fmatrix As; std::vector<giac_double> B_orig,X0; 1005 if (convert(A,As) && convert(b_orig,B_orig) && convert(x0,X0)){ 1006 std::vector<giac_double> res=sparse_conjugate_gradient(As,B_orig,X0,eps,maxiter,contextptr); 1007 return convert(res); 1008 } 1009 } 1010 smatrix As; 1011 if (!convert(A,As)) 1012 return gendimerr(contextptr); 1013 return sparse_conjugate_gradient(As,b_orig,x0,eps,maxiter,contextptr); 1014 } 1015 1016 // Ax=b where A=D+B, Dx_{n+1}=b-B*x_n sparse_jacobi_linsolve(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT)1017 gen sparse_jacobi_linsolve(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT){ 1018 int n=int(A.m.size()); 1019 smatrix B; 1020 vecteur D(n); 1021 vecteur b=*evalf_double(b_orig,1,contextptr)._VECTptr; 1022 for (int i=0;i<n;++i){ 1023 const vecteur & Ai=*A.m[i]._VECTptr; 1024 const std::vector<int> & posi=A.pos[i]; 1025 vecteur Bi; std::vector<int> posB; 1026 Bi.reserve(Ai.size()); posB.reserve(posi.size()); 1027 for (int j=0;j<int(posi.size());++j){ 1028 if (posi[j]==i){ 1029 D[i]=evalf_double(Ai[j],1,contextptr); 1030 } 1031 else { 1032 posB.push_back(posi[j]); 1033 Bi.push_back(evalf_double(Ai[j],1,contextptr)); 1034 } 1035 } 1036 B.m.push_back(Bi); 1037 B.pos.push_back(posB); 1038 } 1039 vecteur tmp(n),xn(x0),prev(n); 1040 gen bn=l2norm(b,contextptr); 1041 for (int i=0;i<maxiter;++i){ 1042 prev=xn; 1043 sparse_mult(B,xn,tmp); 1044 subvecteur(b,tmp,xn); 1045 iterateur jt=xn.begin(),jtend=xn.end(),dt=D.begin(); 1046 for (;jt!=jtend;++jt){ 1047 *jt=*jt / *dt; 1048 } 1049 gen g=l2norm(xn-prev,contextptr)/bn; 1050 if (is_greater(eps,g,contextptr)) 1051 return xn; 1052 } 1053 *logptr(contextptr) << gettext("Warning! Leaving Jacobi iterative algorithm after maximal number of iterations. Check that your matrix is diagonal dominant.") << '\n'; 1054 return xn; 1055 } 1056 l2norm(const std::vector<giac_double> & v)1057 double l2norm(const std::vector<giac_double> & v){ 1058 double res=0; 1059 std::vector<giac_double>::const_iterator it=v.begin(),itend=v.end(); 1060 for (;it!=itend;++it){ 1061 register double tmp=*it; 1062 res+=tmp*tmp; 1063 } 1064 return std::sqrt(res); 1065 } 1066 subvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b,std::vector<giac_double> & c)1067 void subvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b,std::vector<giac_double> & c){ 1068 if (&b==&c){ 1069 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1070 std::vector<giac_double>::const_iterator at=a.begin(); 1071 for (;ct!=ctend;++ct,++at){ 1072 *ct=*at-*ct; 1073 } 1074 return; 1075 } 1076 if (&a==&c){ 1077 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1078 std::vector<giac_double>::const_iterator bt=b.begin(); 1079 for (;ct!=ctend;++ct,++bt){ 1080 *ct -= *bt; 1081 } 1082 return; 1083 } 1084 c.resize(a.size()); 1085 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1086 std::vector<giac_double>::const_iterator at=a.begin(),bt=b.begin(); 1087 for (;ct!=ctend;++at,++bt,++ct){ 1088 *ct = *at-*bt; 1089 } 1090 } 1091 subvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b)1092 std::vector<giac_double> subvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b){ 1093 std::vector<giac_double> c; 1094 subvecteur(a,b,c); 1095 return c; 1096 } 1097 addvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b,std::vector<giac_double> & c)1098 void addvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b,std::vector<giac_double> & c){ 1099 if (&b==&c){ 1100 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1101 std::vector<giac_double>::const_iterator at=a.begin(); 1102 for (;ct!=ctend;++ct,++at){ 1103 *ct=*at+*ct; 1104 } 1105 return; 1106 } 1107 if (&a==&c){ 1108 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1109 std::vector<giac_double>::const_iterator bt=b.begin(); 1110 for (;ct!=ctend;++ct,++bt){ 1111 *ct += *bt; 1112 } 1113 return; 1114 } 1115 c.resize(a.size()); 1116 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1117 std::vector<giac_double>::const_iterator at=a.begin(),bt=b.begin(); 1118 for (;ct!=ctend;++at,++bt,++ct){ 1119 *ct = *at+*bt; 1120 } 1121 } 1122 addvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b)1123 std::vector<giac_double> addvecteur(const std::vector<giac_double> & a,const std::vector<giac_double> & b){ 1124 std::vector<giac_double> c; 1125 addvecteur(a,b,c); 1126 return c; 1127 } 1128 multvecteur(double x,const std::vector<giac_double> & a,std::vector<giac_double> & c)1129 void multvecteur(double x,const std::vector<giac_double> & a,std::vector<giac_double> & c){ 1130 if (&a==&c){ 1131 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1132 for (;ct!=ctend;++ct){ 1133 *ct *= x; 1134 } 1135 return; 1136 } 1137 c.resize(a.size()); 1138 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1139 std::vector<giac_double>::const_iterator at=a.begin(); 1140 for (;ct!=ctend;++at,++ct){ 1141 *ct = x*(*at); 1142 } 1143 } 1144 multvecteur(double x,std::vector<giac_double> & c)1145 void multvecteur(double x,std::vector<giac_double> & c){ 1146 std::vector<giac_double>::iterator ct=c.begin(),ctend=c.end(); 1147 for (;ct!=ctend;++ct){ 1148 *ct *= x; 1149 } 1150 } 1151 multvecteur(double x,const std::vector<giac_double> & b)1152 std::vector<giac_double> multvecteur(double x,const std::vector<giac_double> & b){ 1153 std::vector<giac_double> c(b); 1154 multvecteur(x,c); 1155 return c; 1156 } 1157 1158 // Ax=b where A=D+B, Dx_{n+1}=b-B*x_n 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)1159 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){ 1160 int n=int(A.m.size()); 1161 fmatrix B; 1162 std::vector<giac_double> D(n); 1163 std::vector<giac_double> b=b_orig; 1164 for (int i=0;i<n;++i){ 1165 const std::vector<double> & Ai=A.m[i]; 1166 const std::vector<int> & posi=A.pos[i]; 1167 std::vector<giac_double> Bi; std::vector<int> posB; 1168 Bi.reserve(Ai.size()); posB.reserve(posi.size()); 1169 for (int j=0;j<int(posi.size());++j){ 1170 if (posi[j]==i){ 1171 D[i]=Ai[j]; 1172 } 1173 else { 1174 posB.push_back(posi[j]); 1175 Bi.push_back(Ai[j]); 1176 } 1177 } 1178 B.m.push_back(Bi); 1179 B.pos.push_back(posB); 1180 } 1181 std::vector<giac_double> tmp(n),xn(x0),prev(n); 1182 double bn=l2norm(b); 1183 for (int i=0;i<maxiter;++i){ 1184 prev=xn; 1185 sparse_mult(B,xn,tmp); 1186 subvecteur(b,tmp,xn); 1187 std::vector<giac_double>::iterator jt=xn.begin(),jtend=xn.end(),dt=D.begin(); 1188 for (;jt!=jtend;++jt){ 1189 *jt=*jt / *dt; 1190 } 1191 subvecteur(xn,prev,tmp); 1192 double g=l2norm(tmp)/bn; 1193 if (eps>g) 1194 return xn; 1195 } 1196 *logptr(contextptr) << gettext("Warning! Leaving Jacobi iterative algorithm after maximal number of iterations. Check that your matrix is diagonal dominant.") << '\n'; 1197 return xn; 1198 } 1199 sparse_jacobi_linsolve(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT)1200 gen sparse_jacobi_linsolve(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double eps,int maxiter,GIAC_CONTEXT){ 1201 fmatrix Asf; 1202 std::vector<giac_double> B_orig,X0; 1203 if (convert(A,Asf) && convert(b_orig,B_orig) && convert(x0,X0)){ 1204 std::vector<giac_double> res=sparse_jacobi_linsolve(Asf,B_orig,X0,eps,maxiter,contextptr); 1205 return convert(res); 1206 } 1207 smatrix As; 1208 if (!convert(A,As)) 1209 return gendimerr(contextptr); 1210 return sparse_jacobi_linsolve(As,b_orig,x0,eps,maxiter,contextptr); 1211 } 1212 1213 // Ax=b where A=L+D+U, (D+L)x_{n+1}=b-U*x_n (Gauss-Seidel for omega==1) 1214 // or (L+D/omega)*x_{n+1}=b-(U+D*(1-1/omega))*x_n 1215 // or (-omega*E+D)*x_{n+1}=omega*b+(omega*F+D*(1-omega))*x_n 1216 // or (I-omega*D^-1*E)x_{n+1}=omega*D^-1*b+((1-omega)*I+omega*D^-1*F)*x_n sparse_gauss_seidel_linsolve(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double omega,double eps,int maxiter,GIAC_CONTEXT)1217 gen sparse_gauss_seidel_linsolve(const smatrix & A,const vecteur & b_orig,const vecteur & x0,double omega,double eps,int maxiter,GIAC_CONTEXT){ 1218 int n=int(A.m.size()); 1219 double invomega=1/omega; 1220 smatrix L,U; 1221 vecteur b=*evalf_double(b_orig,1,contextptr)._VECTptr; 1222 for (int i=0;i<n;++i){ 1223 const vecteur & Ai=*A.m[i]._VECTptr; 1224 const std::vector<int> & posi=A.pos[i]; 1225 vecteur Li,Ui; std::vector<int> posL,posU; 1226 for (int j=0;j<int(posi.size());++j){ 1227 if (posi[j]==i) { 1228 posL.push_back(posi[j]); 1229 Li.push_back(invomega*evalf_double(Ai[j],1,contextptr)); 1230 if (invomega!=1){ 1231 posU.push_back(posi[j]); 1232 Ui.push_back((1-invomega)*evalf_double(Ai[j],1,contextptr)); 1233 } 1234 continue; 1235 } 1236 if (posi[j]<i){ 1237 posL.push_back(posi[j]); 1238 Li.push_back(evalf_double(Ai[j],1,contextptr)); 1239 continue; 1240 } 1241 posU.push_back(posi[j]); 1242 Ui.push_back(evalf_double(Ai[j],1,contextptr)); 1243 } 1244 L.m.push_back(Li); 1245 L.pos.push_back(posL); 1246 U.m.push_back(Ui); 1247 U.pos.push_back(posU); 1248 } 1249 vecteur tmp(n),xn(x0),prev(n); 1250 gen bn=l2norm(b,contextptr); 1251 gen_map Lfixme; convert(L,Lfixme); 1252 for (int i=0;i<maxiter;++i){ 1253 prev=xn; 1254 sparse_mult(U,xn,tmp); 1255 subvecteur(b,tmp,tmp); 1256 if (!sparse_linsolve_l(Lfixme,tmp,xn)) 1257 return gensizeerr(contextptr); 1258 gen g=l2norm(xn-prev,contextptr)/bn; 1259 if (is_greater(eps,g,contextptr)) 1260 return xn; 1261 } 1262 *logptr(contextptr) << gettext("Warning! Leaving Gauss-Seidel iterative algorithm after maximal number of iterations. Check that your matrix is diagonal dominant.") << '\n'; 1263 return xn; 1264 } 1265 1266 // Ax=b where A=L+D+U, (D+L)x_{n+1}=b-U*x_n (Gauss-Seidel for omega==1) 1267 // or (L+D/omega)*x_{n+1}=b-(U+D*(1-1/omega))*x_n 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)1268 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){ 1269 int n=int(A.m.size()); 1270 double bn=l2norm(b_orig); 1271 #if 1 1272 // Wikipedia notations 1273 std::vector<giac_double> tmp(n),phiknew(n),phik(x0); 1274 for (int iter=0;iter<maxiter;++iter){ 1275 for (int i=0;i<n;++i){ 1276 giac_double sigma=0,aii=0; 1277 // const std::vector<giac_double> Ai=A.m[i]; 1278 const std::vector<int> & posi=A.pos[i]; 1279 bool ok=false; 1280 std::vector<int>::const_iterator post=posi.begin(),postend=posi.end(); 1281 std::vector<giac_double>::const_iterator at=A.m[i].begin(); 1282 for (;post!=postend;++at,++post){ 1283 int posij=*post; 1284 if (posij==i) 1285 aii=*at; 1286 else 1287 sigma += *at*(posij<i?phiknew[posij]:phik[posij]); 1288 } 1289 if (aii==0) 1290 return std::vector<giac_double>(0); 1291 phiknew[i]=(1-omega)*phik[i]+omega/aii*(b_orig[i]-sigma); 1292 } 1293 subvecteur(phik,phiknew,tmp); 1294 if (l2norm(tmp)<eps*bn){ 1295 if (debug_infolevel) 1296 *logptr(contextptr) << "Convergence criterium reached after " << iter << " iterations, omega=" << omega << '\n'; 1297 return phiknew; 1298 } 1299 phik.swap(phiknew); 1300 } 1301 *logptr(contextptr) << gettext("Warning! Leaving Gauss-Seidel iterative algorithm after maximal number of iterations. Check that your matrix is symetric definite.") << '\n'; 1302 return phik; 1303 #else 1304 std::vector<giac_double> b=b_orig; 1305 double invomega=1/omega; 1306 fmatrix L,U; 1307 for (int i=0;i<n;++i){ 1308 const std::vector<giac_double> Ai=A.m[i]; 1309 const std::vector<int> & posi=A.pos[i]; 1310 std::vector<giac_double> Li,Ui; std::vector<int> posL,posU; 1311 for (int j=0;j<int(posi.size());++j){ 1312 if (posi[j]==i) { 1313 posL.push_back(posi[j]); 1314 Li.push_back(invomega*Ai[j]); 1315 if (invomega!=1){ 1316 posU.push_back(posi[j]); 1317 Ui.push_back((1-invomega)*Ai[j]); 1318 } 1319 continue; 1320 } 1321 if (posi[j]<i){ 1322 posL.push_back(posi[j]); 1323 Li.push_back(Ai[j]); 1324 continue; 1325 } 1326 posU.push_back(posi[j]); 1327 Ui.push_back(Ai[j]); 1328 } 1329 L.m.push_back(Li); 1330 L.pos.push_back(posL); 1331 U.m.push_back(Ui); 1332 U.pos.push_back(posU); 1333 } 1334 std::vector<giac_double> tmp(n),xn(x0),prev(n); 1335 for (int i=0;i<maxiter;++i){ 1336 prev=xn; 1337 sparse_mult(U,xn,tmp); 1338 subvecteur(b,tmp,tmp); 1339 if (!sparse_linsolve_l(L,tmp,xn)) 1340 return std::vector<giac_double>(0); 1341 // CERR << xn << '\n'; 1342 subvecteur(xn,prev,tmp); 1343 double g=l2norm(tmp)/bn; 1344 if (eps>=g){ 1345 if (debug_infolevel) 1346 *logptr(contextptr) << "Convergence criterium reached after " << i << " iterations, omega=" << omega << '\n'; 1347 return xn; 1348 } 1349 } 1350 *logptr(contextptr) << gettext("Warning! Leaving Gauss-Seidel iterative algorithm after maximal number of iterations. Check that your matrix is symetric definite.") << '\n'; 1351 return xn; 1352 #endif 1353 } 1354 sparse_gauss_seidel_linsolve(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double omega,double eps,int maxiter,GIAC_CONTEXT)1355 gen sparse_gauss_seidel_linsolve(const gen_map & A,const vecteur & b_orig,const vecteur & x0,double omega,double eps,int maxiter,GIAC_CONTEXT){ 1356 fmatrix Asf; 1357 std::vector<giac_double> B_orig,X0; 1358 if (convert(A,Asf) && convert(b_orig,B_orig) && convert(x0,X0)){ 1359 std::vector<giac_double> res=sparse_gauss_seidel_linsolve(Asf,B_orig,X0,omega,eps,maxiter,contextptr); 1360 return convert(res); 1361 } 1362 smatrix As; 1363 if (!convert(A,As)) 1364 return gendimerr(contextptr); 1365 return sparse_gauss_seidel_linsolve(As,b_orig,x0,omega,eps,maxiter,contextptr); 1366 } 1367 1368 #ifndef NO_NAMESPACE_GIAC 1369 } 1370 #endif // ndef NO_NAMESPACE_GIAC 1371