1 /**************************************************************************** 2 * VCGLib o o * 3 * Visual and Computer Graphics Library o o * 4 * _ O _ * 5 * Copyright(C) 2004-2016 \/)\/ * 6 * Visual Computing Lab /\/| * 7 * ISTI - Italian National Research Council | * 8 * \ * 9 * All rights reserved. * 10 * * 11 * This program is free software; you can redistribute it and/or modify * 12 * it under the terms of the GNU General Public License as published by * 13 * the Free Software Foundation; either version 2 of the License, or * 14 * (at your option) any later version. * 15 * * 16 * This program is distributed in the hope that it will be useful, * 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * 20 * for more details. * 21 * * 22 ****************************************************************************/ 23 /*************************************************************************** 24 $Log: not supported by cvs2svn $ 25 Revision 1.9 2006/09/11 16:11:39 marfr960 26 Added const to declarations of the overloaded (operators *). 27 Otherwise the * operator would always attempt to convert any type of data passed as an argument to Point3<TYPE> 28 29 Revision 1.8 2006/08/23 15:24:45 marfr960 30 Copy constructor : faster memcpy instead of slow 'for' cycle 31 empty constructor 32 33 Revision 1.7 2006/04/29 10:26:04 fiorin 34 Added some utility methods (swapping of columns and rows, matrix-vector multiplication) 35 36 Revision 1.6 2006/04/11 08:09:35 zifnab1974 37 changes necessary for gcc 3.4.5 on linux 64bit. Please take note of case-sensitivity of filenames 38 39 Revision 1.5 2005/12/12 11:25:00 ganovelli 40 added diagonal matrix, outer produce and namespace 41 42 ***************************************************************************/ 43 44 #ifndef MATRIX_VCGLIB 45 #define MATRIX_VCGLIB 46 47 #include <stdio.h> 48 #include <math.h> 49 #include <memory.h> 50 #include <assert.h> 51 #include <algorithm> 52 #include <vcg/space/point.h> 53 54 namespace vcg{ 55 namespace ndim{ 56 57 /** \addtogroup math */ 58 /* @{ */ 59 60 /*! 61 * This class represent a diagonal <I>m</I>�<I>m</I> matrix. 62 */ 63 64 class MatrixDiagBase{public: 65 virtual const int & Dimension()const =0; 66 virtual float operator[](const int & i)const = 0; 67 }; 68 template<int N, class S> 69 class MatrixDiag: public Point<N,S>, public MatrixDiagBase{ 70 public: Dimension()71 const int & Dimension() const {return N;} MatrixDiag(const Point<N,S> & p)72 MatrixDiag(const Point<N,S>&p):Point<N,S>(p){} 73 }; 74 75 /*! 76 * This class represent a generic <I>m</I>�<I>n</I> matrix. The class is templated over the scalar type field. 77 * @param TYPE (Templete Parameter) Specifies the ScalarType field. 78 */ 79 template<class TYPE> 80 class Matrix 81 { 82 83 public: 84 typedef TYPE ScalarType; 85 86 /*! 87 * Default constructor 88 * All the elements are initialized to zero. 89 * \param m the number of matrix rows 90 * \param n the number of matrix columns 91 */ Matrix(unsigned int m,unsigned int n)92 Matrix(unsigned int m, unsigned int n) 93 { 94 _rows = m; 95 _columns = n; 96 _data = new ScalarType[m*n]; 97 memset(_data, 0, m*n*sizeof(ScalarType)); 98 }; 99 100 /*! 101 * Constructor 102 * The matrix elements are initialized with the values of the elements in \i values. 103 * \param m the number of matrix rows 104 * \param n the number of matrix columns 105 * \param values the values of the matrix elements 106 */ Matrix(unsigned int m,unsigned int n,TYPE * values)107 Matrix(unsigned int m, unsigned int n, TYPE *values) 108 { 109 _rows = m; 110 _columns = n; 111 unsigned int dim = m*n; 112 _data = new ScalarType[dim]; 113 memcpy(_data, values, dim*sizeof(ScalarType)); 114 //unsigned int i; 115 //for (i=0; i<_rows*_columns; i++) 116 // _data[i] = values[i]; 117 }; 118 119 /*! 120 * Empty constructor 121 * Just create the object 122 */ Matrix()123 Matrix() 124 { 125 _rows = 0; 126 _columns = 0; 127 _data = NULL; 128 }; 129 130 /*! 131 * Copy constructor 132 * The matrix elements are initialized with the value of the corresponding element in \i m 133 * \param m the matrix to be copied 134 */ Matrix(const Matrix<TYPE> & m)135 Matrix(const Matrix<TYPE> &m) 136 { 137 _rows = m._rows; 138 _columns = m._columns; 139 _data = new ScalarType[_rows*_columns]; 140 141 unsigned int dim = _rows * _columns; 142 memcpy(_data, m._data, dim * sizeof(ScalarType)); 143 144 // for (unsigned int i=0; i<_rows*_columns; i++) 145 // _data[i] = m._data[i]; 146 }; 147 148 /*! 149 * Default destructor 150 */ ~Matrix()151 ~Matrix() 152 { 153 delete []_data; 154 }; 155 156 /*! 157 * Number of columns 158 */ ColumnsNumber()159 inline unsigned int ColumnsNumber() const 160 { 161 return _columns; 162 }; 163 164 165 /*! 166 * Number of rows 167 */ RowsNumber()168 inline unsigned int RowsNumber() const 169 { 170 return _rows; 171 }; 172 173 /*! 174 * Equality operator. 175 * \param m 176 * \return true iff the matrices have same size and its elements have same values. 177 */ 178 bool operator==(const Matrix<TYPE> &m) const 179 { 180 if (_rows==m._rows && _columns==m._columns) 181 { 182 bool result = true; 183 for (unsigned int i=0; i<_rows*_columns && result; i++) 184 result = (_data[i]==m._data[i]); 185 return result; 186 } 187 return false; 188 }; 189 190 /*! 191 * Inequality operator 192 * \param m 193 * \return true iff the matrices have different size or if their elements have different values. 194 */ 195 bool operator!=(const Matrix<TYPE> &m) const 196 { 197 if (_rows==m._rows && _columns==m._columns) 198 { 199 bool result = false; 200 for (unsigned int i=0; i<_rows*_columns && !result; i++) 201 result = (_data[i]!=m._data[i]); 202 return result; 203 } 204 return true; 205 }; 206 207 /*! 208 * Return the element stored in the <I>i</I>-th rows at the <I>j</I>-th column 209 * \param i the row index 210 * \param j the column index 211 * \return the element 212 */ ElementAt(unsigned int i,unsigned int j)213 inline TYPE ElementAt(unsigned int i, unsigned int j) 214 { 215 assert(i>=0 && i<_rows); 216 assert(j>=0 && j<_columns); 217 return _data[i*_columns+j]; 218 }; 219 220 /*! 221 * Calculate and return the matrix determinant (Laplace) 222 * \return the matrix determinant 223 */ Determinant()224 TYPE Determinant() const 225 { 226 assert(_rows == _columns); 227 switch (_rows) 228 { 229 case 2: 230 { 231 return _data[0]*_data[3]-_data[1]*_data[2]; 232 break; 233 }; 234 case 3: 235 { 236 return _data[0]*(_data[4]*_data[8]-_data[5]*_data[7]) - 237 _data[1]*(_data[3]*_data[8]-_data[5]*_data[6]) + 238 _data[2]*(_data[3]*_data[7]-_data[4]*_data[6]) ; 239 break; 240 }; 241 default: 242 { 243 // da migliorare: si puo' cercare la riga/colonna con maggior numero di zeri 244 ScalarType det = 0; 245 for (unsigned int j=0; j<_columns; j++) 246 if (_data[j]!=0) 247 det += _data[j]*this->Cofactor(0, j); 248 249 return det; 250 } 251 }; 252 }; 253 254 /*! 255 * Return the cofactor <I>A<SUB>i,j</SUB></I> of the <I>a<SUB>i,j</SUB></I> element 256 * \return ... 257 */ Cofactor(unsigned int i,unsigned int j)258 TYPE Cofactor(unsigned int i, unsigned int j) const 259 { 260 assert(_rows == _columns); 261 assert(_rows>2); 262 TYPE *values = new TYPE[(_rows-1)*(_columns-1)]; 263 unsigned int u, v, p, q, s, t; 264 for (u=0, p=0, s=0, t=0; u<_rows; u++, t+=_rows) 265 { 266 if (i==u) 267 continue; 268 269 for (v=0, q=0; v<_columns; v++) 270 { 271 if (j==v) 272 continue; 273 values[s+q] = _data[t+v]; 274 q++; 275 } 276 p++; 277 s+=(_rows-1); 278 } 279 Matrix<TYPE> temp(_rows-1, _columns-1, values); 280 return (pow(TYPE(-1.0), TYPE(i+j))*temp.Determinant()); 281 }; 282 283 /*! 284 * Subscript operator: 285 * \param i the index of the row 286 * \return a reference to the <I>i</I>-th matrix row 287 */ 288 inline TYPE* operator[](const unsigned int i) 289 { 290 assert(i<_rows); 291 return _data + i*_columns; 292 }; 293 294 /*! 295 * Const subscript operator 296 * \param i the index of the row 297 * \return a reference to the <I>i</I>-th matrix row 298 */ 299 inline const TYPE* operator[](const unsigned int i) const 300 { 301 assert(i<_rows); 302 return _data + i*_columns; 303 }; 304 305 /*! 306 * Get the <I>j</I>-th column on the matrix. 307 * \param j the column index. 308 * \return the reference to the column elements. This pointer must be deallocated by the caller. 309 */ GetColumn(const unsigned int j)310 TYPE* GetColumn(const unsigned int j) 311 { 312 assert(j>=0 && j<_columns); 313 ScalarType *v = new ScalarType[_columns]; 314 unsigned int i, p; 315 for (i=0, p=j; i<_rows; i++, p+=_columns) 316 v[i] = _data[p]; 317 return v; 318 }; 319 320 /*! 321 * Get the <I>i</I>-th row on the matrix. 322 * \param i the column index. 323 * \return the reference to the row elements. This pointer must be deallocated by the caller. 324 */ GetRow(const unsigned int i)325 TYPE* GetRow(const unsigned int i) 326 { 327 assert(i>=0 && i<_rows); 328 ScalarType *v = new ScalarType[_rows]; 329 unsigned int j, p; 330 for (j=0, p=i*_columns; j<_columns; j++, p++) 331 v[j] = _data[p]; 332 return v; 333 }; 334 335 /*! 336 * Swaps the values of the elements between the <I>i</I>-th and the <I>j</I>-th column. 337 * \param i the index of the first column 338 * \param j the index of the second column 339 */ SwapColumns(const unsigned int i,const unsigned int j)340 void SwapColumns(const unsigned int i, const unsigned int j) 341 { 342 assert(0<=i && i<_columns); 343 assert(0<=j && j<_columns); 344 if (i==j) 345 return; 346 347 unsigned int r, e0, e1; 348 for (r=0, e0=i, e1=j; r<_rows; r++, e0+=_columns, e1+=_columns) 349 std::swap(_data[e0], _data[e1]); 350 }; 351 352 /*! 353 * Swaps the values of the elements between the <I>i</I>-th and the <I>j</I>-th row. 354 * \param i the index of the first row 355 * \param j the index of the second row 356 */ SwapRows(const unsigned int i,const unsigned int j)357 void SwapRows(const unsigned int i, const unsigned int j) 358 { 359 assert(0<=i && i<_rows); 360 assert(0<=j && j<_rows); 361 if (i==j) 362 return; 363 364 unsigned int r, e0, e1; 365 for (r=0, e0=i*_columns, e1=j*_columns; r<_columns; r++, e0++, e1++) 366 std::swap(_data[e0], _data[e1]); 367 }; 368 369 /*! 370 * Assignment operator 371 * \param m ... 372 */ 373 Matrix<TYPE>& operator=(const Matrix<TYPE> &m) 374 { 375 if (this != &m) 376 { 377 assert(_rows == m._rows); 378 assert(_columns == m._columns); 379 for (unsigned int i=0; i<_rows*_columns; i++) 380 _data[i] = m._data[i]; 381 } 382 return *this; 383 }; 384 385 /*! 386 * Adds a matrix <I>m</I> to this matrix. 387 * \param m reference to matrix to add to this 388 * \return the matrix sum. 389 */ 390 Matrix<TYPE>& operator+=(const Matrix<TYPE> &m) 391 { 392 assert(_rows == m._rows); 393 assert(_columns == m._columns); 394 for (unsigned int i=0; i<_rows*_columns; i++) 395 _data[i] += m._data[i]; 396 return *this; 397 }; 398 399 /*! 400 * Subtracts a matrix <I>m</I> to this matrix. 401 * \param m reference to matrix to subtract 402 * \return the matrix difference. 403 */ 404 Matrix<TYPE>& operator-=(const Matrix<TYPE> &m) 405 { 406 assert(_rows == m._rows); 407 assert(_columns == m._columns); 408 for (unsigned int i=0; i<_rows*_columns; i++) 409 _data[i] -= m._data[i]; 410 return *this; 411 }; 412 413 /*! 414 * (Modifier) Add to each element of this matrix the scalar constant <I>k</I>. 415 * \param k the scalar constant 416 * \return the modified matrix 417 */ 418 Matrix<TYPE>& operator+=(const TYPE k) 419 { 420 for (unsigned int i=0; i<_rows*_columns; i++) 421 _data[i] += k; 422 return *this; 423 }; 424 425 /*! 426 * (Modifier) Subtract from each element of this matrix the scalar constant <I>k</I>. 427 * \param k the scalar constant 428 * \return the modified matrix 429 */ 430 Matrix<TYPE>& operator-=(const TYPE k) 431 { 432 for (unsigned int i=0; i<_rows*_columns; i++) 433 _data[i] -= k; 434 return *this; 435 }; 436 437 /*! 438 * (Modifier) Multiplies each element of this matrix by the scalar constant <I>k</I>. 439 * \param k the scalar constant 440 * \return the modified matrix 441 */ 442 Matrix<TYPE>& operator*=(const TYPE k) 443 { 444 for (unsigned int i=0; i<_rows*_columns; i++) 445 _data[i] *= k; 446 return *this; 447 }; 448 449 /*! 450 * (Modifier) Divides each element of this matrix by the scalar constant <I>k</I>. 451 * \param k the scalar constant 452 * \return the modified matrix 453 */ 454 Matrix<TYPE>& operator/=(const TYPE k) 455 { 456 assert(k!=0); 457 for (unsigned int i=0; i<_rows*_columns; i++) 458 _data[i] /= k; 459 return *this; 460 }; 461 462 /*! 463 * Matrix multiplication: calculates the cross product. 464 * \param m reference to the matrix to multiply by 465 * \return the matrix product 466 */ 467 Matrix<TYPE> operator*(const Matrix<TYPE> &m) const 468 { 469 assert(_columns == m._rows); 470 Matrix<TYPE> result(_rows, m._columns); 471 unsigned int i, j, k, p, q, r; 472 for (i=0, p=0, r=0; i<result._rows; i++, p+=_columns, r+=result._columns) 473 for (j=0; j<result._columns; j++) 474 { 475 ScalarType temp = 0; 476 for (k=0, q=0; k<_columns; k++, q+=m._columns) 477 temp+=(_data[p+k]*m._data[q+j]); 478 result._data[r+j] = temp; 479 } 480 481 return result; 482 }; 483 484 /*! 485 * Matrix-Vector product. Computes the product of the matrix by the vector v. 486 * \param v reference to the vector to multiply by 487 * \return the matrix-vector product. This pointer must be deallocated by the caller 488 */ 489 ScalarType* operator*(const ScalarType v[]) const 490 { 491 ScalarType *result = new ScalarType[_rows]; 492 memset(result, 0, _rows*sizeof(ScalarType)); 493 unsigned int r, c, i; 494 for (r=0, i=0; r<_rows; r++) 495 for (c=0; c<_columns; c++, i++) 496 result[r] += _data[i]*v[c]; 497 498 return result; 499 }; 500 501 /*! 502 * Matrix multiplication: calculates the cross product. 503 * \param reference to the matrix to multiply by 504 * \return the matrix product 505 */ 506 template <int N,int M> DotProduct(Point<N,TYPE> & m,Point<M,TYPE> & result)507 void DotProduct(Point<N,TYPE> &m,Point<M,TYPE> &result) 508 { 509 unsigned int i, j, p, r; 510 for (i=0, p=0, r=0; i<M; i++) 511 { result[i]=0; 512 for (j=0; j<N; j++) 513 result[i]+=(*this)[i][j]*m[j]; 514 } 515 }; 516 517 /*! 518 * Matrix multiplication by a diagonal matrix 519 */ 520 Matrix<TYPE> operator*(const MatrixDiagBase &m) const 521 { 522 assert(_columns == _rows); 523 assert(_columns == m.Dimension()); 524 int i,j; 525 Matrix<TYPE> result(_rows, _columns); 526 527 for (i=0; i<result._rows; i++) 528 for (j=0; j<result._columns; j++) 529 result[i][j]*= m[j]; 530 531 return result; 532 }; 533 /*! 534 * Matrix from outer product. 535 */ 536 template <int N, int M> OuterProduct(const Point<N,TYPE> a,const Point<M,TYPE> b)537 void OuterProduct(const Point<N,TYPE> a, const Point< M,TYPE> b) 538 { 539 assert(N == _rows); 540 assert(M == _columns); 541 Matrix<TYPE> result(_rows,_columns); 542 unsigned int i, j; 543 544 for (i=0; i<result._rows; i++) 545 for (j=0; j<result._columns; j++) 546 (*this)[i][j] = a[i] * b[j]; 547 }; 548 549 550 /*! 551 * Matrix-vector multiplication. 552 * \param reference to the 3-dimensional vector to multiply by 553 * \return the resulting vector 554 */ 555 556 Point3<TYPE> operator*(Point3<TYPE> &p) const 557 { 558 assert(_columns==3 && _rows==3); 559 vcg::Point3<TYPE> result; 560 result[0] = _data[0]*p[0]+_data[1]*p[1]+_data[2]*p[2]; 561 result[1] = _data[3]*p[0]+_data[4]*p[1]+_data[5]*p[2]; 562 result[2] = _data[6]*p[0]+_data[7]*p[1]+_data[8]*p[2]; 563 return result; 564 }; 565 566 567 /*! 568 * Scalar sum. 569 * \param k 570 * \return the resultant matrix 571 */ 572 Matrix<TYPE> operator+(const TYPE k) 573 { 574 Matrix<TYPE> result(_rows, _columns); 575 for (unsigned int i=0; i<_rows*_columns; i++) 576 result._data[i] = _data[i]+k; 577 return result; 578 }; 579 580 /*! 581 * Scalar difference. 582 * \param k 583 * \return the resultant matrix 584 */ 585 Matrix<TYPE> operator-(const TYPE k) 586 { 587 Matrix<TYPE> result(_rows, _columns); 588 for (unsigned int i=0; i<_rows*_columns; i++) 589 result._data[i] = _data[i]-k; 590 return result; 591 }; 592 593 /*! 594 * Negate all matrix elements 595 * \return the modified matrix 596 */ 597 Matrix<TYPE> operator-() const 598 { 599 Matrix<TYPE> result(_rows, _columns, _data); 600 for (unsigned int i=0; i<_columns*_rows; i++) 601 result._data[i] = -1*_data[i]; 602 return result; 603 }; 604 605 /*! 606 * Scalar multiplication. 607 * \param k value to multiply every member by 608 * \return the resultant matrix 609 */ 610 Matrix<TYPE> operator*(const TYPE k) const 611 { 612 Matrix<TYPE> result(_rows, _columns); 613 for (unsigned int i=0; i<_rows*_columns; i++) 614 result._data[i] = _data[i]*k; 615 return result; 616 }; 617 618 /*! 619 * Scalar division. 620 * \param k value to divide every member by 621 * \return the resultant matrix 622 */ 623 Matrix<TYPE> operator/(const TYPE k) 624 { 625 Matrix<TYPE> result(_rows, _columns); 626 for (unsigned int i=0; i<_rows*_columns; i++) 627 result._data[i] = _data[i]/k; 628 return result; 629 }; 630 631 632 /*! 633 * Set all the matrix elements to zero. 634 */ SetZero()635 void SetZero() 636 { 637 for (unsigned int i=0; i<_rows*_columns; i++) 638 _data[i] = ScalarType(0.0); 639 }; 640 641 /*! 642 * Set the matrix to identity. 643 */ SetIdentity()644 void SetIdentity() 645 { 646 assert(_rows==_columns); 647 for (unsigned int i=0; i<_rows; i++) 648 for (unsigned int j=0; j<_columns; j++) 649 _data[i] = (i==j) ? ScalarType(1.0) : ScalarType(0.0f); 650 }; 651 652 /*! 653 * Set the values of <I>j</I>-th column to v[j] 654 * \param j the column index 655 * \param v ... 656 */ SetColumn(const unsigned int j,TYPE * v)657 void SetColumn(const unsigned int j, TYPE* v) 658 { 659 assert(j>=0 && j<_columns); 660 unsigned int i, p; 661 for (i=0, p=j; i<_rows; i++, p+=_columns) 662 _data[p] = v[i]; 663 }; 664 665 /*! 666 * Set the elements of the <I>i</I>-th row to v[j] 667 * \param i the row index 668 * \param v ... 669 */ SetRow(const unsigned int i,TYPE * v)670 void SetRow(const unsigned int i, TYPE* v) 671 { 672 assert(i>=0 && i<_rows); 673 unsigned int j, p; 674 for (j=0, p=i*_rows; j<_columns; j++, p++) 675 _data[p] = v[j]; 676 }; 677 678 /*! 679 * Set the diagonal elements <I>v<SUB>i,i</SUB></I> to v[i] 680 * \param v 681 */ SetDiagonal(TYPE * v)682 void SetDiagonal(TYPE *v) 683 { 684 assert(_rows == _columns); 685 for (unsigned int i=0, p=0; i<_rows; i++, p+=_rows) 686 _data[p+i] = v[i]; 687 }; 688 689 /*! 690 * Resize the current matrix. 691 * \param m the number of matrix rows. 692 * \param n the number of matrix columns. 693 */ Resize(const unsigned int m,const unsigned int n)694 void Resize(const unsigned int m, const unsigned int n) 695 { 696 assert(m>=2); 697 assert(n>=2); 698 _rows = m; 699 _columns = n; 700 delete []_data; 701 _data = new ScalarType[m*n]; 702 for (unsigned int i=0; i<m*n; i++) 703 _data[i] = 0; 704 }; 705 706 707 /*! 708 * Matrix transposition operation: set the current matrix to its transpose 709 */ Transpose()710 void Transpose() 711 { 712 ScalarType *temp = new ScalarType[_rows*_columns]; 713 unsigned int i, j, p, q; 714 for (i=0, p=0; i<_rows; i++, p+=_columns) 715 for (j=0, q=0; j<_columns; j++, q+=_rows) 716 temp[q+i] = _data[p+j]; 717 718 std::swap(_columns, _rows); 719 std::swap(_data, temp); 720 delete []temp; 721 }; 722 723 // for the transistion to eigen transpose()724 Matrix transpose() 725 { 726 Matrix res = *this; 727 res.Transpose(); 728 return res; 729 } transposeInPlace()730 void transposeInPlace() { Transpose(); } 731 // for the transistion to eigen 732 733 /*! 734 * Print all matrix elements 735 */ Dump()736 void Dump() 737 { 738 unsigned int i, j, p; 739 for (i=0, p=0; i<_rows; i++, p+=_columns) 740 { 741 printf("[\t"); 742 for (j=0; j<_columns; j++) 743 printf("%f\t", _data[p+j]); 744 745 printf("]\n"); 746 } 747 printf("\n"); 748 }; 749 750 protected: 751 /// the number of matrix rows 752 unsigned int _rows; 753 754 /// the number of matrix rows 755 unsigned int _columns; 756 757 /// the matrix elements 758 ScalarType *_data; 759 }; 760 761 typedef vcg::ndim::Matrix<double> MatrixMNd; 762 typedef vcg::ndim::Matrix<float> MatrixMNf; 763 764 /*! @} */ 765 766 // template <class MatrixType> 767 // void Invert(MatrixType & m){ 768 // typedef typename MatrixType::ScalarType X; 769 // X *diag; 770 // diag = new X [m.ColumnsNumber()]; 771 772 // MatrixType res(m.RowsNumber(),m.ColumnsNumber()); 773 // vcg::SingularValueDecomposition<MatrixType > (m,&diag[0],res,LeaveUnsorted,50 ); 774 // m.Transpose(); 775 // // prodotto per la diagonale 776 // unsigned int i,j; 777 // for (i=0; i<m.RowsNumber(); i++) 778 // for (j=0; j<m.ColumnsNumber(); j++) 779 // res[i][j]/= diag[j]; 780 // m = res *m; 781 // } 782 783 } 784 }; // end of namespace 785 786 #endif 787