1/* linbox/matrix/matrix-domain.inl 2 * Copyright (C) 2002 Bradford Hovinen 3 * 4 * Written by Bradford Hovinen <bghovinen@math.uwaterloo.ca> 5 * 6 * ------------------------------------ 7 * 8 * 9 * ========LICENCE======== 10 * This file is part of the library LinBox. 11 * 12 * LinBox is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU Lesser General Public 14 * License as published by the Free Software Foundation; either 15 * version 2.1 of the License, or (at your option) any later version. 16 * 17 * This library is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 20 * Lesser General Public License for more details. 21 * 22 * You should have received a copy of the GNU Lesser General Public 23 * License along with this library; if not, write to the Free Software 24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 25 * ========LICENCE======== 26 *. 27 */ 28 29#ifndef __LINBOX_matrix_domain_INL 30#define __LINBOX_matrix_domain_INL 31 32#include <linbox/linbox-config.h> 33#include "linbox/matrix/transpose-matrix.h" 34#include "linbox/blackbox/dif.h" 35 36namespace LinBox 37{ 38 39 template<class Field> 40 template <class Matrix1, class Matrix2> 41 Matrix1 &MatrixDomain<Field>::copyRow (Matrix1 &A, const Matrix2 &B) const 42 { 43 linbox_check (A.rowdim () == B.rowdim ()); 44 linbox_check (A.coldim () == B.coldim ()); 45 46 typename Matrix1::RowIterator i; 47 typename Matrix2::ConstRowIterator j; 48 49 i = A.rowBegin (); 50 j = B.rowBegin (); 51 52 for (; i != A.rowEnd (); ++i, ++j) 53 _VD.copy (*i, *j); 54 55 return A; 56 } 57 58 template<class Field> 59 template <class Matrix1, class Matrix2> 60 Matrix1 &MatrixDomain<Field>::copyCol (Matrix1 &A, const Matrix2 &B) const 61 { 62 linbox_check (A.rowdim () == B.rowdim ()); 63 linbox_check (A.coldim () == B.coldim ()); 64 65 typename Matrix1::ColIterator i; 66 typename Matrix2::ConstColIterator j; 67 68 i = A.colBegin (); 69 j = B.colBegin (); 70 71 for (; i != A.colEnd (); ++i, ++j) 72 _VD.copy (*i, *j); 73 74 return A; 75 } 76 77 template<class Field> 78 template <class Matrix1, class Matrix2> 79 bool MatrixDomain<Field>::areEqualBB (const Matrix1 &A, const Matrix2 &B) const 80 { Dif<Matrix1, Matrix2> C(A, B); return isZero(C); } 81 82 template<class Field> 83 template <class Matrix1, class Matrix2> 84 bool MatrixDomain<Field>::areEqualRow (const Matrix1 &A, const Matrix2 &B) const 85 { 86 linbox_check (A.rowdim () == B.rowdim ()); 87 linbox_check (A.coldim () == B.coldim ()); 88 89 typename Matrix1::ConstRowIterator i; 90 typename Matrix2::ConstRowIterator j; 91 92 i = A.rowBegin (); 93 j = B.rowBegin (); 94 95 for (; i != A.rowEnd (); ++i, ++j) 96 if (!_VD.areEqual (*i, *j)) 97 return false; 98 99 return true; 100 } 101 102 template<class Field> 103 template <class Matrix1, class Matrix2> 104 bool MatrixDomain<Field>::areEqualCol (const Matrix1 &A, const Matrix2 &B) const 105 { 106 linbox_check (A.rowdim () == B.rowdim ()); 107 linbox_check (A.coldim () == B.coldim ()); 108 109 typename Matrix1::ConstColIterator i; 110 typename Matrix2::ConstColIterator j; 111 112 i = A.colBegin (); 113 j = B.colBegin (); 114 115 for (; i != A.colEnd (); ++i, ++j) 116 if (!_VD.areEqual (*i, *j)) 117 return false; 118 119 return true; 120 } 121 122 template<class Field> 123 template <class Matrix_> 124 bool MatrixDomain<Field>::isZeroBB (const Matrix_ &A) const 125 { 126 VectorDomain<Field> VD(A.field()); 127 std::vector<typename Field::Element> x(A.coldim()), y(A.rowdim()); 128 VD.random(x); 129 A.apply(y, x); 130 return VD.isZero(y); 131 } 132 133 template<class Field> 134 template <class Matrix_> 135 bool MatrixDomain<Field>::isZeroRow (const Matrix_ &A) const 136 { 137 typename Matrix_::ConstRowIterator i; 138 139 i = A.rowBegin (); 140 141 for (; i != A.rowEnd (); ++i) 142 if (!_VD.isZero (*i)) 143 return false; 144 145 return true; 146 } 147 148 template<class Field> 149 template <class Matrix_> 150 bool MatrixDomain<Field>::isZeroCol (const Matrix_ &A) const 151 { 152 typename Matrix::ConstColIterator i; 153 154 i = A.colBegin (); 155 156 for (; i != A.colEnd (); ++i) 157 if (!_VD.isZero (*i)) 158 return false; 159 160 return true; 161 } 162 163 template<class Field> 164 template <class Matrix1, class Matrix2, class Matrix3> 165 Matrix1 &MatrixDomain<Field>::addRow (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 166 { 167 linbox_check (A.rowdim () == B.rowdim ()); 168 linbox_check (A.coldim () == B.coldim ()); 169 linbox_check (A.rowdim () == C.rowdim ()); 170 linbox_check (A.coldim () == C.coldim ()); 171 172 typename Matrix1::RowIterator i; 173 typename Matrix2::ConstRowIterator j; 174 typename Matrix3::ConstRowIterator k; 175 176 i = C.rowBegin (); 177 j = A.rowBegin (); 178 k = B.rowBegin (); 179 180 for (; i != C.rowEnd (); ++i, ++j, ++k) 181 _VD.add (*i, *j, *k); 182 183 return C; 184 } 185 186 template<class Field> 187 template <class Matrix1, class Matrix2, class Matrix3> 188 Matrix1 &MatrixDomain<Field>::addCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 189 { 190 linbox_check (A.rowdim () == B.rowdim ()); 191 linbox_check (A.coldim () == B.coldim ()); 192 linbox_check (A.rowdim () == C.rowdim ()); 193 linbox_check (A.coldim () == C.coldim ()); 194 195 typename Matrix1::ColIterator i; 196 typename Matrix2::ConstColIterator j; 197 typename Matrix3::ConstColIterator k; 198 199 i = C.colBegin (); 200 j = A.colBegin (); 201 k = B.colBegin (); 202 203 for (; i != C.colEnd (); ++i, ++j, ++k) 204 _VD.add (*i, *j, *k); 205 206 return C; 207 } 208 209 template<class Field> 210 template <class Matrix1, class Matrix2> 211 Matrix1 &MatrixDomain<Field>::addinRow (Matrix1 &A, const Matrix2 &B) const 212 { 213 linbox_check (A.rowdim () == B.rowdim ()); 214 linbox_check (A.coldim () == B.coldim ()); 215 216 typename Matrix1::RowIterator i; 217 typename Matrix2::ConstRowIterator j; 218 219 i = A.rowBegin (); 220 j = B.rowBegin (); 221 222 for (; i != A.rowEnd (); ++i, ++j) 223 _VD.addin (*i, *j); 224 225 return A; 226 } 227 228 template<class Field> 229 template <class Matrix1, class Matrix2> 230 Matrix1 &MatrixDomain<Field>::addinCol (Matrix1 &A, const Matrix2 &B) const 231 { 232 linbox_check (A.rowdim () == B.rowdim ()); 233 linbox_check (A.coldim () == B.coldim ()); 234 235 typename Matrix1::ColIterator i; 236 typename Matrix2::ConstColIterator j; 237 238 i = A.colBegin (); 239 j = B.colBegin (); 240 241 for (; i != A.colEnd (); ++i, ++j) 242 _VD.addin (*i, *j); 243 244 return A; 245 } 246 247 template<class Field> 248 template <class Matrix1, class Matrix2, class Matrix3> 249 Matrix1 &MatrixDomain<Field>::subRow (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 250 { 251 linbox_check (A.rowdim () == B.rowdim ()); 252 linbox_check (A.coldim () == B.coldim ()); 253 linbox_check (A.rowdim () == C.rowdim ()); 254 linbox_check (A.coldim () == C.coldim ()); 255 256 typename Matrix1::RowIterator i; 257 typename Matrix2::ConstRowIterator j; 258 typename Matrix3::ConstRowIterator k; 259 260 i = C.rowBegin (); 261 j = A.rowBegin (); 262 k = B.rowBegin (); 263 264 for (; i != C.rowEnd (); ++i, ++j, ++k) 265 _VD.sub (*i, *j, *k); 266 267 return C; 268 } 269 270 template<class Field> 271 template <class Matrix1, class Matrix2, class Matrix3> 272 Matrix1 &MatrixDomain<Field>::subCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 273 { 274 linbox_check (A.rowdim () == B.rowdim ()); 275 linbox_check (A.coldim () == B.coldim ()); 276 linbox_check (A.rowdim () == C.rowdim ()); 277 linbox_check (A.coldim () == C.coldim ()); 278 279 typename Matrix1::ColIterator i; 280 typename Matrix2::ConstColIterator j; 281 typename Matrix3::ConstColIterator k; 282 283 i = C.colBegin (); 284 j = A.colBegin (); 285 k = B.colBegin (); 286 287 for (; i != C.colEnd (); ++i, ++j, ++k) 288 _VD.sub (*i, *j, *k); 289 290 return C; 291 } 292 293 template<class Field> 294 template <class Matrix1, class Matrix2> 295 Matrix1 &MatrixDomain<Field>::subinRow (Matrix1 &A, const Matrix2 &B) const 296 { 297 linbox_check (A.rowdim () == B.rowdim ()); 298 linbox_check (A.coldim () == B.coldim ()); 299 300 typename Matrix1::RowIterator i; 301 typename Matrix2::ConstRowIterator j; 302 303 i = A.rowBegin (); 304 j = B.rowBegin (); 305 306 for (; i != A.rowEnd (); ++i, ++j) 307 _VD.subin (*i, *j); 308 309 return A; 310 } 311 312 template<class Field> 313 template <class Matrix1, class Matrix2> 314 Matrix1 &MatrixDomain<Field>::subinCol (Matrix1 &A, const Matrix2 &B) const 315 { 316 linbox_check (A.rowdim () == B.rowdim ()); 317 linbox_check (A.coldim () == B.coldim ()); 318 319 typename Matrix1::ColIterator i; 320 typename Matrix2::ConstColIterator j; 321 322 i = A.colBegin (); 323 j = B.colBegin (); 324 325 for (; i != A.colEnd (); ++i, ++j) 326 _VD.subin (*i, *j); 327 328 return A; 329 } 330 331 template<class Field> 332 template <class Matrix1, class Matrix2> 333 Matrix1 &MatrixDomain<Field>::negRow (Matrix1 &A, const Matrix2 &B) const 334 { 335 linbox_check (A.rowdim () == B.rowdim ()); 336 linbox_check (A.coldim () == B.coldim ()); 337 338 typename Matrix1::RowIterator i; 339 typename Matrix2::ConstRowIterator j; 340 341 i = A.rowBegin (); 342 j = B.rowBegin (); 343 344 for (; i != A.rowEnd (); ++i, ++j) 345 _VD.neg (*i, *j); 346 347 return A; 348 } 349 350 template<class Field> 351 template <class Matrix1, class Matrix2> 352 Matrix1 &MatrixDomain<Field>::negCol (Matrix1 &A, const Matrix2 &B) const 353 { 354 linbox_check (A.rowdim () == B.rowdim ()); 355 linbox_check (A.coldim () == B.coldim ()); 356 357 typename Matrix1::ColIterator i; 358 typename Matrix2::ConstColIterator j; 359 360 i = A.colBegin (); 361 j = B.colBegin (); 362 363 for (; i != A.colEnd (); ++i, ++j) 364 _VD.neg (*i, *j); 365 366 return A; 367 } 368 369 template<class Field> 370 template <class Matrix_> 371 Matrix_ &MatrixDomain<Field>::neginRow (Matrix_ &A) const 372 { 373 typename Matrix_::RowIterator i; 374 375 for (i = A.rowBegin (); i != A.rowEnd (); ++i) 376 _VD.negin (*i); 377 378 return A; 379 } 380 381 template<class Field> 382 template <class Matrix_> 383 Matrix_ &MatrixDomain<Field>::neginCol (Matrix_ &A) const 384 { 385 typename Matrix_::ColIterator i; 386 387 for (i = A.colBegin (); i != A.colEnd (); ++i) 388 _VD.negin (*i); 389 390 return A; 391 } 392 393 template<class Field> 394 template <class Matrix1, class Matrix2, class Matrix3> 395 Matrix1 &MatrixDomain<Field>::mulRowRowCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 396 { 397 linbox_check (A.coldim () == B.rowdim ()); 398 linbox_check (A.rowdim () == C.rowdim ()); 399 linbox_check (B.coldim () == C.coldim ()); 400 401 typename Matrix2::ConstRowIterator i; 402 typename Matrix3::ConstColIterator j; 403 typename Matrix1::RowIterator l1; 404 typename Matrix1::Row::iterator l2; 405 406 for (i = A.rowBegin (), l1 = C.rowBegin (); i != A.rowEnd (); ++i, ++l1) 407 for (j = B.colBegin (), l2 = l1->begin (); j != B.colEnd (); ++j, ++l2) 408 _VD.dot (*l2, *i, *j); 409 410 return C; 411 } 412 413 template<class Field> 414 template <class Matrix1, class Matrix2, class Matrix3> 415 Matrix1 &MatrixDomain<Field>::mulColRowCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 416 { 417 linbox_check (A.coldim () == B.rowdim ()); 418 linbox_check (A.rowdim () == C.rowdim ()); 419 linbox_check (B.coldim () == C.coldim ()); 420 421 typename Matrix2::ConstRowIterator i; 422 typename Matrix3::ConstColIterator j; 423 typename Matrix1::ColIterator l1; 424 typename Matrix1::Col::iterator l2; 425 426 for (j = B.colBegin (), l1 = C.colBegin (); j != B.colEnd (); ++j, ++l1) 427 for (i = A.rowBegin (), l2 = l1->begin (); i != A.rowEnd (); ++i, ++l2) 428 _VD.dot (*l2, *i, *j); 429 430 return C; 431 } 432 433 template<class Field> 434 template <class Matrix1, class Matrix2, class Matrix3> 435 Matrix1 &MatrixDomain<Field>::mulRowRowRow (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 436 { 437 linbox_check (A.coldim () == B.rowdim ()); 438 linbox_check (A.rowdim () == C.rowdim ()); 439 linbox_check (B.coldim () == C.coldim ()); 440 441 typename Matrix2::ConstRowIterator i = A.rowBegin (); 442 typename Matrix1::RowIterator j = C.rowBegin (); 443 444 TransposeMatrix<const Matrix3> BT (B); 445 446 for (; i != A.rowEnd (); ++i, ++j) 447 vectorMul (*j, BT, *i); 448 449 return C; 450 } 451 452 template<class Field> 453 template <class Matrix1, class Matrix2, class Matrix3> 454 Matrix1 &MatrixDomain<Field>::mulColColCol (Matrix1 &C, const Matrix2 &A, const Matrix3 &B) const 455 { 456 linbox_check (A.coldim () == B.rowdim ()); 457 linbox_check (A.rowdim () == C.rowdim ()); 458 linbox_check (B.coldim () == C.coldim ()); 459 460 typename Matrix3::ConstColIterator i = B.colBegin (); 461 typename Matrix1::ColIterator j = C.colBegin (); 462 463 for (; i != B.colEnd (); ++i, ++j) 464 vectorMul (*j, A, *i); 465 466 return C; 467 } 468 469 template<class Field> 470 template <class Matrix1, class Matrix2> 471 Matrix2 &MatrixDomain<Field>::leftMulin (const Matrix1 &A, Matrix2 &B) const 472 { 473 linbox_check (A.rowdim () == A.coldim ()); 474 linbox_check (A.coldim () == B.rowdim ()); 475 476 DenseVector t (field(),A.rowdim ()); 477 478 typename Matrix2::ColIterator i; 479 typename Matrix1::ConstRowIterator j; 480 481 typename DenseVector::iterator k; 482 483 for (i = B.colBegin (); i != B.colEnd (); ++i) { 484 for (j = A.rowBegin (), k = t.begin (); j != A.rowEnd (); ++j, ++k) 485 _VD.dot (*k, *j, *i); 486 487 _VD.copy (*i, t); 488 } 489 490 return B; 491 } 492 493 template<class Field> 494 template <class Matrix1, class Matrix2> 495 Matrix1 &MatrixDomain<Field>::rightMulin (Matrix1 &A, const Matrix2 &B) const 496 { 497 linbox_check (A.coldim () == B.rowdim ()); 498 linbox_check (B.rowdim () == B.coldim ()); 499 500 DenseVector t (field(),B.coldim ()); 501 502 typename Matrix1::RowIterator i; 503 typename Matrix2::ConstColIterator j; 504 505 typename DenseVector::iterator k; 506 507 for (i = A.rowBegin (); i != A.rowEnd (); ++i) { 508 for (j = B.colBegin (), k = t.begin (); j != B.colEnd (); ++j, ++k) 509 _VD.dot (*k, *i, *j); 510 511 _VD.copy (*i, t); 512 } 513 514 return A; 515 } 516 517 template<class Field> 518 template <class Matrix1, class Matrix2> 519 Matrix1 &MatrixDomain<Field>::mulRow (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a) const 520 { 521 linbox_check (C.rowdim () == B.rowdim ()); 522 linbox_check (C.coldim () == B.coldim ()); 523 524 typename Matrix1::RowIterator i; 525 typename Matrix2::ConstRowIterator j; 526 527 i = C.rowBegin (); 528 j = B.rowBegin (); 529 530 for (; i != C.rowEnd (); ++i, ++j) 531 _VD.mul (*i, *j, a); 532 533 return C; 534 } 535 536 template<class Field> 537 template <class Matrix1, class Matrix2> 538 Matrix1 &MatrixDomain<Field>::mulCol (Matrix1 &C, const Matrix2 &B, const typename Field::Element &a) const 539 { 540 linbox_check (C.rowdim () == B.rowdim ()); 541 linbox_check (C.coldim () == B.coldim ()); 542 543 typename Matrix1::ColIterator i; 544 typename Matrix2::ConstColIterator j; 545 546 i = C.colBegin (); 547 j = B.colBegin (); 548 549 for (; i != C.colEnd (); ++i, ++j) 550 _VD.mul (*i, *j, a); 551 552 return C; 553 } 554 555 template<class Field> 556 template <class Matrix_> 557 Matrix_ &MatrixDomain<Field>::mulinRow (Matrix_ &B, const typename Field::Element &a) const 558 { 559 typename Matrix_::RowIterator i; 560 561 for (i = B.rowBegin (); i != B.rowEnd (); ++i) 562 _VD.mulin (*i, a); 563 564 return B; 565 } 566 567 template<class Field> 568 template <class Matrix_> 569 Matrix_ &MatrixDomain<Field>::mulinCol (Matrix_ &B, const typename Field::Element &a) const 570 { 571 typename Matrix_::ColIterator i; 572 573 for (i = B.colBegin (); i != B.colEnd (); ++i) 574 _VD.mulin (*i, a); 575 576 return B; 577 } 578 579 template<class Field> 580 template <class Matrix1, class Matrix2, class Matrix3> 581 Matrix1 &MatrixDomain<Field>::axpyinRowRowCol (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const 582 { 583 linbox_check (A.coldim () == X.rowdim ()); 584 linbox_check (A.rowdim () == Y.rowdim ()); 585 linbox_check (X.coldim () == Y.coldim ()); 586 587 typename Matrix2::ConstRowIterator i; 588 typename Matrix3::ConstColIterator j; 589 typename Matrix1::RowIterator l1; 590 typename Matrix1::Row::iterator l2; 591 592 typename Field::Element t; 593 594 for (i = A.rowBegin (), l1 = Y.rowBegin (); i != A.rowEnd (); ++i, ++l1) { 595 for (j = X.colBegin (), l2 = l1->begin (); j != X.colEnd (); ++j, ++l2) { 596 _VD.dot (t, *i, *j); 597 field().addin (*l2, t); 598 } 599 } 600 601 return Y; 602 } 603 604 template<class Field> 605 template <class Matrix1, class Matrix2, class Matrix3> 606 Matrix1 &MatrixDomain<Field>::axpyinColRowCol (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const 607 { 608 linbox_check (A.coldim () == X.rowdim ()); 609 linbox_check (A.rowdim () == Y.rowdim ()); 610 linbox_check (X.coldim () == Y.coldim ()); 611 612 typename Matrix2::ConstRowIterator i; 613 typename Matrix3::ConstColIterator j; 614 typename Matrix1::ColIterator l1; 615 typename Matrix1::Col::iterator l2; 616 617 typename Field::Element t; 618 619 for (j = X.colBegin (), l1 = Y.colBegin (); j != X.colEnd (); ++j, ++l1) { 620 for (i = A.rowBegin (), l2 = l1->begin (); i != A.rowEnd (); ++i, ++l2) { 621 _VD.dot (t, *i, *j); 622 field().addin (*l2, t); 623 } 624 } 625 626 return Y; 627 } 628 629 template<class Field> 630 template <class Matrix1, class Matrix2, class Matrix3> 631 Matrix1 &MatrixDomain<Field>::axpyinRowRowRow (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const 632 { 633 linbox_check (A.coldim () == X.rowdim ()); 634 linbox_check (A.rowdim () == Y.rowdim ()); 635 linbox_check (X.coldim () == Y.coldim ()); 636 637 DenseVector t (field(),X.coldim ()); 638 639 typename Matrix2::ConstRowIterator i = A.rowBegin (); 640 typename Matrix1::RowIterator j = Y.rowBegin (); 641 642 TransposeMatrix<const Matrix3> XT (X); 643 644 for (; i != A.rowEnd (); ++i, ++j) { 645 vectorMul (t, XT, *i); 646 _VD.addin (*j, t); 647 } 648 649 return Y; 650 } 651 652 template<class Field> 653 template <class Matrix1, class Matrix2, class Matrix3> 654 Matrix1 &MatrixDomain<Field>::axpyinColColCol (Matrix1 &Y, const Matrix2 &A, const Matrix3 &X) const 655 { 656 linbox_check (A.coldim () == X.rowdim ()); 657 linbox_check (A.rowdim () == Y.rowdim ()); 658 linbox_check (X.coldim () == Y.coldim ()); 659 660 DenseVector t (field(),A.rowdim ()); 661 662 typename Matrix3::ConstColIterator i = X.colBegin (); 663 typename Matrix1::ColIterator j = Y.colBegin (); 664 665 for (; i != X.colEnd (); ++i, ++j) { 666 vectorMul (t, A, *i); 667 _VD.addin (*j, t); 668 } 669 670 return Y; 671 } 672 673 template<class Field> 674 template <class Vector1, class Matrix_, class Vector2> 675 Vector1 &MatrixDomain<Field>::mulRowSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 676 VectorCategories::DenseVectorTag) const 677 { 678 linbox_check (A.coldim () == v.size ()); 679 linbox_check (A.rowdim () == w.size ()); 680 681 typename Matrix_::ConstRowIterator i = A.rowBegin (); 682 typename Vector1::iterator j = w.begin (); 683 684 // JGD 02.09.2008 : when sizes differ 685 // A must decide if dot is possible, not w 686 // for (; j != w.end (); ++j, ++i) 687 // _VD.dot (*j, v, *i); 688 for (; i != A.rowEnd (); ++j, ++i) { 689 linbox_check(j != w.end()); 690 this->_VD.dot (*j, v, *i); 691 } 692 693 return w; 694 } 695 696 template<class Field> 697 template <class Vector1, class Matrix_, class Vector2> 698 Vector1 &MatrixDomain<Field>::mulRowSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 699 VectorCategories::SparseSequenceVectorTag) const 700 { 701 typename Matrix_::ConstRowIterator i = A.rowBegin (); 702 typename Field::Element t; 703 unsigned int idx = 0; 704 705 w.clear (); 706 707 for (; i != A.rowEnd (); ++i, ++idx) { 708 _VD.dot (t, v, *i); 709 710 if (!field().isZero (t)) 711 w.push_back (std::pair<size_t, typename Field::Element> (idx, t)); 712 } 713 714 return w; 715 } 716 717 template<class Field> 718 template <class Vector1, class Matrix_, class Vector2> 719 Vector1 &MatrixDomain<Field>::mulRowSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 720 VectorCategories::SparseAssociativeVectorTag) const 721 { 722 typename Matrix_::ConstRowIterator i = A.rowBegin (); 723 typename Field::Element t; 724 unsigned int idx = 0; 725 726 w.clear (); 727 728 for (; i != A.rowEnd (); ++i, ++idx) { 729 _VD.dot (t, v, *i); 730 731 if (!field().isZero (t)) 732 w[idx] = t; 733 } 734 735 return w; 736 } 737 738 template<class Field> 739 template <class Vector1, class Matrix_, class Vector2> 740 Vector1 &MatrixDomain<Field>::mulRowSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 741 VectorCategories::SparseParallelVectorTag) const 742 { 743 typename Matrix_::ConstRowIterator i = A.rowBegin (); 744 typename Field::Element t; 745 unsigned int idx = 0; 746 747 w.first.clear (); 748 w.second.clear (); 749 750 for (; i != A.rowEnd (); ++i, ++idx) { 751 _VD.dot (t, v, *i); 752 753 if (!field().isZero (t)) { 754 w.first.push_back (idx); 755 w.second.push_back (t); 756 } 757 } 758 759 return w; 760 } 761 762 template<class Field> 763 template <class Vector1, class Matrix_, class Vector2> 764 Vector1 &MVProductDomain<Field>::mulColDense 765 (const VectorDomain<Field> &VD, Vector1 &w, const Matrix_ &A, const Vector2 &v) const 766 { 767 linbox_check (A.coldim () == v.size ()); 768 linbox_check (A.rowdim () == w.size ()); 769 770 typename Matrix_::ConstColIterator i = A.colBegin (); 771 typename Vector2::const_iterator j = v.begin (); 772 773 VD.subin (w, w); 774 775 for (; j != v.end (); ++j, ++i) 776 VD.axpyin (w, *j, *i); 777 778 return w; 779 } 780 781 template<class Field> 782 template <class Vector1, class Matrix_, class Vector2> 783 Vector1 &MatrixDomain<Field>::mulColSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 784 VectorCategories::DenseVectorTag, 785 VectorCategories::SparseSequenceVectorTag) const 786 { 787 linbox_check (A.rowdim () == w.size ()); 788 789 typename Vector2::const_iterator j = v.begin (); 790 791 _VD.subin (w, w); 792 793 for (; j != v.end (); ++j) { 794 typename Matrix_::ConstColIterator i = A.colBegin () + j->first; 795 _VD.axpyin (w, j->second, *i); 796 } 797 798 return w; 799 } 800 801 template<class Field> 802 template <class Vector1, class Matrix_, class Vector2> 803 Vector1 &MatrixDomain<Field>::mulColSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 804 VectorCategories::DenseVectorTag, 805 VectorCategories::SparseAssociativeVectorTag) const 806 { 807 linbox_check (A.rowdim () == w.size ()); 808 809 typename Vector2::const_iterator j = v.begin (); 810 811 _VD.subin (w, w); 812 813 for (; j != v.end (); ++j) { 814 typename Matrix_::ConstColIterator i = A.colBegin () + j->first; 815 _VD.axpyin (w, j->second, *i); 816 } 817 818 return w; 819 } 820 821 template<class Field> 822 template <class Vector1, class Matrix_, class Vector2> 823 Vector1 &MatrixDomain<Field>::mulColSpecialized (Vector1 &w, const Matrix_ &A, const Vector2 &v, 824 VectorCategories::DenseVectorTag, 825 VectorCategories::SparseParallelVectorTag) const 826 { 827 linbox_check (A.rowdim () == w.size ()); 828 829 typename Vector2::first_type::const_iterator j_idx = v.first.begin (); 830 typename Vector2::second_type::const_iterator j_elt = v.second.begin (); 831 832 _VD.subin (w, w); 833 834 for (; j_idx != v.first.end (); ++j_idx, ++j_elt) { 835 typename Matrix_::ConstColIterator i = A.colBegin () + (ptrdiff_t)*j_idx; 836 _VD.axpyin (w, *j_elt, *i); 837 } 838 839 return w; 840 } 841 842 template<class Field> 843 template <class Vector1, class Matrix_, class Vector2> 844 Vector1 &MatrixDomain<Field>::axpyinRowSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 845 VectorCategories::DenseVectorTag) const 846 { 847 linbox_check (A.coldim () == x.size ()); 848 linbox_check (A.rowdim () == y.size ()); 849 850 typename Matrix_::ConstRowIterator i = A.rowBegin (); 851 typename Vector1::iterator j = y.begin (); 852 853 typename Field::Element t; 854 855 for (; j != y.end (); ++j, ++i) { 856 _VD.dot (t, x, *i); 857 field().addin (*j, t); 858 } 859 860 return y; 861 } 862 863 template<class Field> 864 template <class Vector1, class Matrix_, class Vector2> 865 Vector1 &MatrixDomain<Field>::axpyinRowSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 866 VectorCategories::SparseSequenceVectorTag) const 867 { 868 DenseVector t1 (field(),A.coldim ()), t2 (field(),A.rowdim ()); 869 870 _VD.copy (t1, x); 871 _VD.copy (t2, y); 872 axpyin (t2, A, t1); 873 _VD.copy (y, t2); 874 875 return y; 876 } 877 878 template<class Field> 879 template <class Vector1, class Matrix_, class Vector2> 880 Vector1 &MatrixDomain<Field>::axpyinRowSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 881 VectorCategories::SparseAssociativeVectorTag) const 882 { 883 DenseVector t1 (field(),A.coldim ()), t2 (field(),A.rowdim ()); 884 885 _VD.copy (t1, x); 886 _VD.copy (t2, y); 887 axpyin (t2, A, t1); 888 _VD.copy (y, t2); 889 890 return y; 891 } 892 893 template<class Field> 894 template <class Vector1, class Matrix_, class Vector2> 895 Vector1 &MatrixDomain<Field>::axpyinRowSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 896 VectorCategories::SparseParallelVectorTag) const 897 { 898 DenseVector t1 (field(),A.coldim ()), t2 (field(),A.rowdim ()); 899 900 _VD.copy (t1, x); 901 _VD.copy (t2, y); 902 axpyin (t2, A, t1); 903 _VD.copy (y, t2); 904 905 return y; 906 } 907 908 template<class Field> 909 template <class Vector1, class Matrix_, class Vector2> 910 Vector1 &MatrixDomain<Field>::axpyinColSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 911 VectorCategories::DenseVectorTag) const 912 { 913 linbox_check (A.coldim () == x.size ()); 914 linbox_check (A.rowdim () == y.size ()); 915 916 typename Matrix_::ConstColIterator i = A.colBegin (); 917 typename Vector2::const_iterator j = x.begin (); 918 919 for (; j != x.end (); ++j, ++i) 920 _VD.axpyin (y, *j, *i); 921 922 return y; 923 } 924 925 template<class Field> 926 template <class Vector1, class Matrix_, class Vector2> 927 Vector1 &MatrixDomain<Field>::axpyinColSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 928 VectorCategories::SparseSequenceVectorTag) const 929 { 930 typename Matrix_::ConstColIterator i = A.colBegin (); 931 typename Vector2::iterator j = x.begin (); 932 933 934 DenseVector t (field(),A.rowdim ()); 935 936 _VD.copy (t, y); 937 938 while (j != x.end ()) { 939 int diff; 940 _VD.axpyin (t, j->second, *i); 941 diff = j->first; ++j; 942 diff -= j->first; 943 i -= diff; 944 } 945 946 _VD.copy (y, t); 947 948 return y; 949 } 950 951 template<class Field> 952 template <class Vector1, class Matrix_, class Vector2> 953 Vector1 &MatrixDomain<Field>::axpyinColSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 954 VectorCategories::SparseAssociativeVectorTag) const 955 { 956 typename Matrix_::ConstColIterator i = A.colBegin (); 957 typename Vector2::iterator j = x.begin (); 958 959 960 DenseVector t (field(),A.rowdim ()); 961 962 _VD.copy (t, y); 963 964 while (j != x.end ()) { 965 int diff; 966 _VD.axpyin (t, j->second, *i); 967 diff = j->first; ++j; 968 diff -= j->first; 969 i -= diff; 970 } 971 972 _VD.copy (y, t); 973 974 return y; 975 } 976 977 template<class Field> 978 template <class Vector1, class Matrix_, class Vector2> 979 Vector1 &MatrixDomain<Field>::axpyinColSpecialized (Vector1 &y, const Matrix_ &A, const Vector2 &x, 980 VectorCategories::SparseParallelVectorTag) const 981 { 982 typename Matrix_::ConstColIterator i = A.colBegin (); 983 typename Vector2::iterator j_idx = x.first.begin (); 984 typename Vector2::iterator j_elt = x.second.begin (); 985 986 987 DenseVector t (field(),A.rowdim ()); 988 989 _VD.copy (t, y); 990 991 for (; j_idx != x.first.end (); ++j_elt) { 992 int diff; 993 _VD.axpyin (t, *j_elt, *i); 994 diff = *j_idx; ++j_idx; 995 diff -= *j_idx; 996 i -= diff; 997 } 998 999 _VD.copy (y, t); 1000 1001 return y; 1002 } 1003 1004 template<class Field> 1005 template <class Matrix1, class Blackbox, class Matrix2> 1006 Matrix1 &MatrixDomain<Field>::blackboxMulLeft (Matrix1 &C, const Blackbox &A, const Matrix2 &B) const 1007 { 1008 linbox_check (A.coldim () == B.rowdim ()); 1009 linbox_check (A.rowdim () == C.rowdim ()); 1010 linbox_check (B.coldim () == C.coldim ()); 1011 1012 typename Matrix1::ColIterator i = C.colBegin (); 1013 typename Matrix2::ConstColIterator j = B.colBegin (); 1014 1015 for (; i != C.colEnd (); ++i, ++j) 1016 A.apply (*i, *j); 1017 1018 return C; 1019 } 1020 1021 template<class Field> 1022 template <class Matrix1, class Matrix2, class Blackbox> 1023 Matrix1 &MatrixDomain<Field>::blackboxMulRight (Matrix1 &C, const Matrix2 &A, const Blackbox &B) const 1024 { 1025 linbox_check (A.coldim () == B.rowdim ()); 1026 linbox_check (A.rowdim () == C.rowdim ()); 1027 linbox_check (B.coldim () == C.coldim ()); 1028 1029 typename Matrix1::RowIterator i = C.rowBegin (); 1030 typename Matrix2::ConstRowIterator j = A.rowBegin (); 1031 1032 for (; i != C.rowEnd (); ++i, ++j) 1033 B.applyTranspose (*i, *j); 1034 1035 return C; 1036 } 1037 1038 template<class Field> 1039 template <class Matrix_, class Iterator> 1040 Matrix_ &MatrixDomain<Field>::permuteRowsByRow (Matrix_ &A, 1041 Iterator P_start, 1042 Iterator P_end) const 1043 { 1044 Iterator i; 1045 typename Matrix_::RowIterator j, k; 1046 1047 for (i = P_start; i != P_end; ++i) { 1048 j = A.rowBegin () + i->first; 1049 k = A.rowBegin () + i->second; 1050 1051 _VD.swap (*j, *k); 1052 } 1053 1054 return A; 1055 } 1056 1057 template<class Field> 1058 template <class Matrix_, class Iterator> 1059 Matrix_ &MatrixDomain<Field>::permuteRowsByCol (Matrix_ &A, 1060 Iterator P_start, 1061 Iterator P_end) const 1062 { 1063 typename Matrix_::ColIterator j; 1064 1065 for (j = A.colBegin (); j != A.colEnd (); ++j) 1066 _VD.permute (*j, P_start, P_end); 1067 1068 return A; 1069 } 1070 1071 template<class Field> 1072 template <class Matrix_, class Iterator> 1073 Matrix_ &MatrixDomain<Field>::permuteColsByRow (Matrix_ &A, 1074 Iterator P_start, 1075 Iterator P_end) const 1076 { 1077 typename Matrix_::RowIterator j; 1078 1079 for (j = A.rowBegin (); j != A.rowEnd (); ++j) 1080 _VD.permute (*j, P_start, P_end); 1081 1082 return A; 1083 } 1084 1085 template<class Field> 1086 template <class Matrix_, class Iterator> 1087 Matrix_ &MatrixDomain<Field>::permuteColsByCol (Matrix_ &A, 1088 Iterator P_start, 1089 Iterator P_end) const 1090 { 1091 Iterator i; 1092 typename Matrix_::ColIterator j, k; 1093 1094 for (i = P_start; i != P_end; ++i) { 1095 j = A.colBegin () + i->first; 1096 k = A.colBegin () + i->second; 1097 1098 _VD.swap (*j, *k); 1099 } 1100 1101 return A; 1102 } 1103 1104 /* FIXME: These methods are undocumented, and I'm unclear what they are supposed 1105 * to do 1106 */ 1107 1108#if 0 1109 1110 /*M1<-M2**k; 1111 */ 1112 template<class Matrix1, class Matrix2> 1113 Matrix1& MatrixDomain::pow_apply(Matrix1& M1, const Matrix2& M2, unsigned long int k) const 1114 { 1115 linbox_check((M1.rowdim()==M1.coldim())&& 1116 (M2.rowdim()==M2.coldim())&& 1117 (M1.rowdim()==M2.rowdim())); 1118 1119 1120 typename Matrix1::Iterator p=M1.Begin(); 1121 for(;p!=M1.End();++p) 1122 M1.field().assign(*p,M1.field().zero); 1123 for(p=M1.Begin();p<M1.End();) 1124 { 1125 M1.field().assign(*p,M1.field().one); 1126 p=p+M1.rowdim()+1; 1127 } 1128 1129 1130 for(int i=0;i<k;++i) 1131 mulin_R(M1,M2); 1132 1133 return M1; 1134 } 1135 1136 1137 template<class Matrix1, class Matrix2> 1138 Matrix1& MatrixDomain::pow_horn(Matrix1& M1, const Matrix2& M2, unsigned long int k) const 1139 { 1140 linbox_check((M1.rowdim()==M1.coldim())&& 1141 (M2.rowdim()==M2.coldim())&& 1142 (M1.rowdim()==M2.rowdim())); 1143 1144 if(k==0) 1145 { 1146 typename Matrix1::Iterator p=M1.Begin(); 1147 for(;p!=M1.End();++p) 1148 M1.field().assign(*p,M1.field().zero); 1149 for(p=M1.Begin();p<M1.End();) 1150 { 1151 M1.field().assign(*p,M1.field().one); 1152 p+=M1.rowdim()+1; 1153 } 1154 return M1; 1155 } 1156 1157 typename Matrix1::Iterator p1; 1158 typename Matrix2::ConstIterator p2; 1159 for(p1=M1.Begin(),p2=M2.Begin();p1!=M1.End();++p1,++p2) 1160 M1.field().assign(*p1,*p2); 1161 1162 std::vector<bool> bit; 1163 bit.reserve(sizeof(unsigned long)*4); 1164 while(k>0) 1165 { 1166 bit.push_back(k%2); 1167 k/=2; 1168 }; 1169 1170 1171 std::vector<bool>::reverse_iterator p=bit.rbegin(); 1172 ++p; 1173 Matrix1 temp(M1); 1174 for(;p!=bit.rend();++p) 1175 { 1176 temp=M1; 1177 mulin_L(M1,temp); 1178 if(*p) 1179 mulin_L(M1,M2); 1180 1181 } 1182 1183 return M1; 1184 } 1185 1186#endif 1187 1188} // namespace LinBox 1189 1190#endif // __LINBOX_matrix_domain_INL 1191 1192// Local Variables: 1193// mode: C++ 1194// tab-width: 4 1195// indent-tabs-mode: nil 1196// c-basic-offset: 4 1197// End: 1198// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 1199