1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "InterpKernelMatrixTools.hxx"
22 #include "InterpKernelException.hxx"
23 #include "InterpKernelAutoPtr.hxx"
28 namespace INTERP_KERNEL
31 * Computes the dot product of two vectors.
32 * This routine uses unrolled loops for increments equal to one.
36 * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
37 * LINPACK User's Guide,
39 * ISBN13: 978-0-898711-72-1,
42 * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
43 * Basic Linear Algebra Subprograms for Fortran Usage,
45 * ACM Transactions on Mathematical Software,
46 * Volume 5, Number 3, September 1979, pages 308-323.
48 * \param [in] n the number of entries in the vectors.
49 * \param [in] dx the first vector.
50 * \param [in] incx the increment between successive entries in \a dx.
51 * \param [in] dy the second vector.
52 * \param [in] incy the increment between successive entries in \a dy.
53 * \return the sum of the product of the corresponding entries of \a dx and \a dy.
55 double ddot(int n, const double *dx, int incx, const double *dy, int incy)
61 // Code for unequal increments or equal increments not equal to 1.
62 if(incx!=1 || incy!=1)
73 for(i=0;i<n;i++,ix+=incx,iy+=incy)
76 //Code for both increments equal to 1.
83 dtemp+=dx[i]*dy[i]+dx[i+1]*dy[i+1]+dx[i+2]*dy[i+2]+dx[i+3]*dy[i+3]+dx[i+4]*dy[i+4];
89 void dscal(int n, double sa, double *x, int incx)
116 for(i=0;i<n;i++,ix+=incx)
121 void daxpy(int n, double da, const double *dx, int incx, double *dy, int incy)
128 // Code for unequal increments or equal increments not equal to 1.
129 if(incx!=1 || incy!=1)
141 for(i=0;i<n;i++,ix+=incx,iy+=incy)
142 dy[iy]=dy[iy]+da*dx[ix];
144 // Code for both increments equal to 1.
149 dy[i]=dy[i]+da*dx[i];
152 dy[i ]=dy[i ]+da*dx[i ];
153 dy[i+1]=dy[i+1]+da*dx[i+1];
154 dy[i+2]=dy[i+2]+da*dx[i+2];
155 dy[i+3]=dy[i+3]+da*dx[i+3];
160 double r8_abs(double x)
168 void dswap(int n, double *x, int incx, double *y, int incy)
174 else if(incx==1 && incy==1)
178 { temp=x[i]; x[i]=y[i]; y[i]=temp; }
181 temp=x[i]; x[i]=y[i]; y[i]=temp;
182 temp=x[i+1]; x[i+1]=y[i+1]; y[i+1]=temp;
183 temp=x[i+2]; x[i+2]=y[i+2]; y[i+2]=temp;
196 for(i=0;i<n;i++,ix+=incx,iy+=incy)
198 temp=x[ix]; x[ix]=y[iy];
204 * Finds the index of the vector element of maximum absolute value.
205 * \warning This index is a 1-based index, not a 0-based index !
208 * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
209 * LINPACK User's Guide,
211 * ISBN13: 978-0-898711-72-1,
214 * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
215 * Basic Linear Algebra Subprograms for Fortran Usage,
217 * ACM Transactions on Mathematical Software,
218 * Volume 5, Number 3, September 1979, pages 308-323.
220 * \param [in] n the number of entries in the vector.
221 * \param [in] dx the vector to be examined.
222 * \param [in] incx the increment between successive entries of SX.
223 * \return the index of the element of maximum absolute value (in C convention).
225 int idamax(int n, const double *dx, int incx)
230 if ( n < 1 || incx <= 0 )
240 if(dmax<r8_abs(dx[i]))
254 if(dmax<r8_abs(dx[ix]))
269 * factors a real general matrix.
273 * Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
274 * LINPACK User's Guide,
275 * SIAM, (Society for Industrial and Applied Mathematics),
276 * 3600 University City Science Center,
277 * Philadelphia, PA, 19104-2688.
282 * \param [in,out] a a matrix of size LDA*N. On intput, the matrix to be factored.
283 * On output, an upper triangular matrix and the multipliers used to obtain
284 * it. The factorization can be written A=L*U, where L is a product of
285 * permutation and unit lower triangular matrices, and U is upper triangular.
287 * \param [in] lda the leading dimension of \a a.
288 * \param [in] n the order of the matrix \a a.
289 * \param [out] ipvt the pivot indices of size \a n.
290 * \return the singularity indicator.
292 * - K, if U(K-1,K-1) == 0. This is not an error condition for this subroutine,
293 * but it does indicate that DGESL or DGEDI will divide by zero if called.
295 int dgefa(double *a, int lda, int n, int *ipvt)
300 // Gaussian elimination with partial pivoting.
301 for(int k=0;k<n-1;k++)
303 // Find L=pivot index.
304 l=idamax(n-k-1,a+k+k*lda,1)+k;
306 // Zero pivot implies this column already triangularized.
312 //Interchange if necessary.
316 a[l+k*lda]=a[k+k*lda];
319 // Compute multipliers.
321 dscal(n-k-1,t,a+k+1+k*lda,1);
322 // Row elimination with column indexing.
323 for(int j=k+1;j<n;j++)
328 a[l+j*lda]=a[k+j*lda];
331 daxpy(n-k-1,t,a+k+1+k*lda,1,a+k+1+j*lda,1);
335 if(a[n-1+(n-1)*lda]==0.0)
341 * - Purpose: computes inverse (job=1) of a matrix factored by dgefa.
343 * - Discussion: A division by zero will occur if the input factor contains
344 * a zero on the diagonal. It will not occur if the subroutines are called correctly
345 * and dgfa has set info to 0.
349 * Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
350 * LINPACK User's Guide,
351 * SIAM, (Society for Industrial and Applied Mathematics),
352 * 3600 University City Science Center,
353 * Philadelphia, PA, 19104-2688.
358 * \param [in,out] a matrix of size lda*n, on input, the LU factor information, as output of dgfa. On output, the inverse matrix.
359 * \param [in] lda, the leading dimension of the array \a a.
360 * \param [in] n, the order of the matrix \a a.
361 * \param [in] ipvt, the pivot vector from dgfa.
362 * \param [in,out] work a work array of size at least equal to \a n.
364 void dgedi(double *a, int lda, int n, const int *ipvt, double *work)
369 a[k+k*lda]=1.0/a[k+k*lda];
371 dscal(k,t,a+0+k*lda,1);
372 for(int j=k+1;j<n;j++)
376 daxpy(k+1,t,a+0+k*lda,1,a+0+j*lda,1);
379 // Form inverse(U) * inverse(L).
380 for(int k=n-2;k>=0;k--)
382 for(int i=k+1;i<n;i++)
388 for(int j=k+1;j<n;j++)
391 daxpy(n,t,a+0+j*lda,1,a+0+k*lda,1);
395 dswap(n,a+0+k*lda,1,a+0+l*lda,1);
399 void matrixProduct(const double *A, int n1, int p1, const double *B, int n2, int p2, double *C)
403 std::ostringstream oss; oss << "matrixProduct : the size of input matrix are not coherent the nb of cols of input matrix #0 is " << p1 << " whereas the number of rows of input matrix #1 is " << n2 << " !";
404 throw INTERP_KERNEL::Exception(oss.str().c_str());
406 for(int i=0;i<n1;i++)
408 for(int j=0;j<p2;j++)
411 for(int k=0;k<p1;k++)
412 C[i*p2+j]+=A[i*p1+k]*B[k*p2+j];
417 void inverseMatrix(const double *A, int n, double *iA)
419 INTERP_KERNEL::AutoPtr<int> ipvt=new int[n];
420 INTERP_KERNEL::AutoPtr<double> work=new double[n*n];
421 std::copy(A,A+n*n,iA);
423 dgedi(iA,n,n,ipvt,work);