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 "InterpKernelAutoPtr.hxx"
26 namespace INTERP_KERNEL
29 * Computes the dot product of two vectors.
30 * This routine uses unrolled loops for increments equal to one.
34 * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
35 * LINPACK User's Guide,
37 * ISBN13: 978-0-898711-72-1,
40 * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
41 * Basic Linear Algebra Subprograms for Fortran Usage,
43 * ACM Transactions on Mathematical Software,
44 * Volume 5, Number 3, September 1979, pages 308-323.
46 * \param [in] n the number of entries in the vectors.
47 * \param [in] dx the first vector.
48 * \param [in] incx the increment between successive entries in \a dx.
49 * \param [in] dy the second vector.
50 * \param [in] incy the increment between successive entries in \a dy.
51 * \return the sum of the product of the corresponding entries of \a dx and \a dy.
53 double ddot(int n, const double *dx, int incx, const double *dy, int incy)
59 // Code for unequal increments or equal increments not equal to 1.
60 if(incx!=1 || incy!=1)
71 for(i=0;i<n;i++,ix+=incx,iy+=incy)
74 //Code for both increments equal to 1.
81 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];
87 void dscal(int n, double sa, double *x, int incx)
114 for(i=0;i<n;i++,ix+=incx)
119 void daxpy(int n, double da, const double *dx, int incx, double *dy, int incy)
126 // Code for unequal increments or equal increments not equal to 1.
127 if(incx!=1 || incy!=1)
139 for(i=0;i<n;i++,ix+=incx,iy+=incy)
140 dy[iy]=dy[iy]+da*dx[ix];
142 // Code for both increments equal to 1.
147 dy[i]=dy[i]+da*dx[i];
150 dy[i ]=dy[i ]+da*dx[i ];
151 dy[i+1]=dy[i+1]+da*dx[i+1];
152 dy[i+2]=dy[i+2]+da*dx[i+2];
153 dy[i+3]=dy[i+3]+da*dx[i+3];
158 double r8_abs(double x)
166 void dswap(int n, double *x, int incx, double *y, int incy)
172 else if(incx==1 && incy==1)
176 { temp=x[i]; x[i]=y[i]; y[i]=temp; }
179 temp=x[i]; x[i]=y[i]; y[i]=temp;
180 temp=x[i+1]; x[i+1]=y[i+1]; y[i+1]=temp;
181 temp=x[i+2]; x[i+2]=y[i+2]; y[i+2]=temp;
194 for(i=0;i<n;i++,ix+=incx,iy+=incy)
196 temp=x[ix]; x[ix]=y[iy];
202 * Finds the index of the vector element of maximum absolute value.
203 * \warning This index is a 1-based index, not a 0-based index !
206 * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
207 * LINPACK User's Guide,
209 * ISBN13: 978-0-898711-72-1,
212 * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
213 * Basic Linear Algebra Subprograms for Fortran Usage,
215 * ACM Transactions on Mathematical Software,
216 * Volume 5, Number 3, September 1979, pages 308-323.
218 * \param [in] n the number of entries in the vector.
219 * \param [in] dx the vector to be examined.
220 * \param [in] incx the increment between successive entries of SX.
221 * \return the index of the element of maximum absolute value (in C convention).
223 int idamax(int n, const double *dx, int incx)
228 if ( n < 1 || incx <= 0 )
238 if(dmax<r8_abs(dx[i]))
252 if(dmax<r8_abs(dx[ix]))
267 * factors a real general matrix.
271 * Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
272 * LINPACK User's Guide,
273 * SIAM, (Society for Industrial and Applied Mathematics),
274 * 3600 University City Science Center,
275 * Philadelphia, PA, 19104-2688.
280 * \param [in,out] a a matrix of size LDA*N. On intput, the matrix to be factored.
281 * On output, an upper triangular matrix and the multipliers used to obtain
282 * it. The factorization can be written A=L*U, where L is a product of
283 * permutation and unit lower triangular matrices, and U is upper triangular.
285 * \param [in] lda the leading dimension of \a a.
286 * \param [in] n the order of the matrix \a a.
287 * \param [out] ipvt the pivot indices of size \a n.
288 * \return the singularity indicator.
290 * - K, if U(K-1,K-1) == 0. This is not an error condition for this subroutine,
291 * but it does indicate that DGESL or DGEDI will divide by zero if called.
293 int dgefa(double *a, int lda, int n, int *ipvt)
298 // Gaussian elimination with partial pivoting.
299 for(int k=0;k<n-1;k++)
301 // Find L=pivot index.
302 l=idamax(n-k-1,a+k+k*lda,1)+k;
304 // Zero pivot implies this column already triangularized.
310 //Interchange if necessary.
314 a[l+k*lda]=a[k+k*lda];
317 // Compute multipliers.
319 dscal(n-k-1,t,a+k+1+k*lda,1);
320 // Row elimination with column indexing.
321 for(int j=k+1;j<n;j++)
326 a[l+j*lda]=a[k+j*lda];
329 daxpy(n-k-1,t,a+k+1+k*lda,1,a+k+1+j*lda,1);
333 if(a[n-1+(n-1)*lda]==0.0)
339 * - Purpose: computes inverse (job=1) of a matrix factored by dgefa.
341 * - Discussion: A division by zero will occur if the input factor contains
342 * a zero on the diagonal. It will not occur if the subroutines are called correctly
343 * and dgfa has set info to 0.
347 * Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
348 * LINPACK User's Guide,
349 * SIAM, (Society for Industrial and Applied Mathematics),
350 * 3600 University City Science Center,
351 * Philadelphia, PA, 19104-2688.
356 * \param [in,out] a matrix of size lda*n, on input, the LU factor information, as output of dgfa. On output, the inverse matrix.
357 * \param [in] lda, the leading dimension of the array \a a.
358 * \param [in] n, the order of the matrix \a a.
359 * \param [in] ipvt, the pivot vector from dgfa.
360 * \param [in,out] work a work array of size at least equal to \a n.
362 void dgedi(double *a, int lda, int n, const int *ipvt, double *work)
367 a[k+k*lda]=1.0/a[k+k*lda];
369 dscal(k,t,a+0+k*lda,1);
370 for(int j=k+1;j<n;j++)
374 daxpy(k+1,t,a+0+k*lda,1,a+0+j*lda,1);
377 // Form inverse(U) * inverse(L).
378 for(int k=n-2;k>=0;k--)
380 for(int i=k+1;i<n;i++)
386 for(int j=k+1;j<n;j++)
389 daxpy(n,t,a+0+j*lda,1,a+0+k*lda,1);
393 dswap(n,a+0+k*lda,1,a+0+l*lda,1);
397 void matrixProduct(const double *A, int n1, int p1, const double *B, int n2, int p2, double *C)
399 for(int i=0;i<n1;i++)
401 for(int j=0;j<p2;j++)
404 for(int k=0;k<p1;k++)
405 C[i*p2+j]+=A[i*p1+k]*B[k*p2+j];
410 void inverseMatrix(const double *A, int n, double *iA)
412 INTERP_KERNEL::AutoPtr<int> ipvt=new int[n];
413 INTERP_KERNEL::AutoPtr<double> work=new double[n*n];
414 std::copy(A,A+n*n,iA);
416 dgedi(iA,n,n,ipvt,work);