1 // Copyright (C) 2007-2012 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
20 #include "InterpKernelMatrixTools.hxx"
21 #include "InterpKernelAutoPtr.hxx"
25 namespace INTERP_KERNEL
28 * Computes the dot product of two vectors.
29 * This routine uses unrolled loops for increments equal to one.
33 * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
34 * LINPACK User's Guide,
36 * ISBN13: 978-0-898711-72-1,
39 * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
40 * Basic Linear Algebra Subprograms for Fortran Usage,
42 * ACM Transactions on Mathematical Software,
43 * Volume 5, Number 3, September 1979, pages 308-323.
45 * \param [in] n the number of entries in the vectors.
46 * \param [in] dx the first vector.
47 * \param [in] incx the increment between successive entries in \a dx.
48 * \param [in] dy the second vector.
49 * \param [in] incy the increment between successive entries in \a dy.
50 * \return the sum of the product of the corresponding entries of \a dx and \a dy.
52 double ddot(int n, const double *dx, int incx, const double *dy, int incy)
58 // Code for unequal increments or equal increments not equal to 1.
59 if(incx!=1 || incy!=1)
70 for(i=0;i<n;i++,ix+=incx,iy+=incy)
73 //Code for both increments equal to 1.
80 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];
86 void dscal(int n, double sa, double *x, int incx)
113 for(i=0;i<n;i++,ix+=incx)
118 void daxpy(int n, double da, const double *dx, int incx, double *dy, int incy)
125 // Code for unequal increments or equal increments not equal to 1.
126 if(incx!=1 || incy!=1)
138 for(i=0;i<n;i++,ix+=incx,iy+=incy)
139 dy[iy]=dy[iy]+da*dx[ix];
141 // Code for both increments equal to 1.
146 dy[i]=dy[i]+da*dx[i];
149 dy[i ]=dy[i ]+da*dx[i ];
150 dy[i+1]=dy[i+1]+da*dx[i+1];
151 dy[i+2]=dy[i+2]+da*dx[i+2];
152 dy[i+3]=dy[i+3]+da*dx[i+3];
157 double r8_abs(double x)
165 void dswap(int n, double *x, int incx, double *y, int incy)
171 else if(incx==1 && incy==1)
175 { temp=x[i]; x[i]=y[i]; y[i]=temp; }
178 temp=x[i]; x[i]=y[i]; y[i]=temp;
179 temp=x[i+1]; x[i+1]=y[i+1]; y[i+1]=temp;
180 temp=x[i+2]; x[i+2]=y[i+2]; y[i+2]=temp;
193 for(i=0;i<n;i++,ix+=incx,iy+=incy)
195 temp=x[ix]; x[ix]=y[iy];
201 * Finds the index of the vector element of maximum absolute value.
202 * \warning This index is a 1-based index, not a 0-based index !
205 * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
206 * LINPACK User's Guide,
208 * ISBN13: 978-0-898711-72-1,
211 * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
212 * Basic Linear Algebra Subprograms for Fortran Usage,
214 * ACM Transactions on Mathematical Software,
215 * Volume 5, Number 3, September 1979, pages 308-323.
217 * \param [in] n the number of entries in the vector.
218 * \param [in] dx the vector to be examined.
219 * \param [in] incx the increment between successive entries of SX.
220 * \return the index of the element of maximum absolute value (in C convention).
222 int idamax(int n, const double *dx, int incx)
227 if ( n < 1 || incx <= 0 )
237 if(dmax<r8_abs(dx[i]))
251 if(dmax<r8_abs(dx[ix]))
266 * factors a real general matrix.
270 * Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
271 * LINPACK User's Guide,
272 * SIAM, (Society for Industrial and Applied Mathematics),
273 * 3600 University City Science Center,
274 * Philadelphia, PA, 19104-2688.
279 * \param [in,out] a a matrix of size LDA*N. On intput, the matrix to be factored.
280 * On output, an upper triangular matrix and the multipliers used to obtain
281 * it. The factorization can be written A=L*U, where L is a product of
282 * permutation and unit lower triangular matrices, and U is upper triangular.
284 * \param [in] lda the leading dimension of \a a.
285 * \param [in] n the order of the matrix \a a.
286 * \param [out] ipvt the pivot indices of size \a n.
287 * \return the singularity indicator.
289 * - K, if U(K-1,K-1) == 0. This is not an error condition for this subroutine,
290 * but it does indicate that DGESL or DGEDI will divide by zero if called.
292 int dgefa(double *a, int lda, int n, int *ipvt)
297 // Gaussian elimination with partial pivoting.
298 for(int k=0;k<n-1;k++)
300 // Find L=pivot index.
301 l=idamax(n-k-1,a+k+k*lda,1)+k;
303 // Zero pivot implies this column already triangularized.
309 //Interchange if necessary.
313 a[l+k*lda]=a[k+k*lda];
316 // Compute multipliers.
318 dscal(n-k-1,t,a+k+1+k*lda,1);
319 // Row elimination with column indexing.
320 for(int j=k+1;j<n;j++)
325 a[l+j*lda]=a[k+j*lda];
328 daxpy(n-k-1,t,a+k+1+k*lda,1,a+k+1+j*lda,1);
332 if(a[n-1+(n-1)*lda]==0.0)
338 * - Purpose: computes inverse (job=1) of a matrix factored by dgefa.
340 * - Discussion: A division by zero will occur if the input factor contains
341 * a zero on the diagonal. It will not occur if the subroutines are called correctly
342 * and dgfa has set info to 0.
346 * Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
347 * LINPACK User's Guide,
348 * SIAM, (Society for Industrial and Applied Mathematics),
349 * 3600 University City Science Center,
350 * Philadelphia, PA, 19104-2688.
355 * \param [in,out] a matrix of size lda*n, on input, the LU factor information, as output of dgfa. On output, the inverse matrix.
356 * \param [in] lda, the leading dimension of the array \a a.
357 * \param [in] n, the order of the matrix \a a.
358 * \param [in] ipvt, the pivot vector from dgfa.
359 * \param [in,out] work a work array of size at least equal to \a n.
361 void dgedi(double *a, int lda, int n, const int *ipvt, double *work)
366 a[k+k*lda]=1.0/a[k+k*lda];
368 dscal(k,t,a+0+k*lda,1);
369 for(int j=k+1;j<n;j++)
373 daxpy(k+1,t,a+0+k*lda,1,a+0+j*lda,1);
376 // Form inverse(U) * inverse(L).
377 for(int k=n-2;k>=0;k--)
379 for(int i=k+1;i<n;i++)
385 for(int j=k+1;j<n;j++)
388 daxpy(n,t,a+0+j*lda,1,a+0+k*lda,1);
392 dswap(n,a+0+k*lda,1,a+0+l*lda,1);
396 void matrixProduct(const double *A, int n1, int p1, const double *B, int n2, int p2, double *C)
398 for(int i=0;i<n1;i++)
400 for(int j=0;j<p2;j++)
403 for(int k=0;k<p1;k++)
404 C[i*p2+j]+=A[i*p1+k]*B[k*p2+j];
409 void inverseMatrix(const double *A, int n, double *iA)
411 INTERP_KERNEL::AutoPtr<int> ipvt=new int[n];
412 INTERP_KERNEL::AutoPtr<double> work=new double[n*n];
413 std::copy(A,A+n*n,iA);
415 dgedi(iA,n,n,ipvt,work);