2 #include "gmx_lapack.h"
3 #include "lapack_limits.h"
6 F77_FUNC(dgetrf,DGETRF)(int *m,
16 double minusone = -1.0;
24 mindim = (*m < *n) ? *m : *n;
26 if(DGETRF_BLOCKSIZE>=mindim) {
29 F77_FUNC(dgetf2,DGETF2)(m,n,a,lda,ipiv,info);
35 for(j=1;j<=mindim;j+=DGETRF_BLOCKSIZE) {
36 jb = ( DGETRF_BLOCKSIZE < (mindim-j+1)) ? DGETRF_BLOCKSIZE : (mindim-j+1);
37 /* factor diag. and subdiag blocks and test for singularity */
39 F77_FUNC(dgetf2,DGETF2)(&k,&jb,&(a[(j-1)*(*lda)+(j-1)]),lda,&(ipiv[j-1]),&iinfo);
41 if(*info==0 && iinfo>0)
42 *info = iinfo + j - 1;
44 /* adjust pivot indices */
45 k = (*m < (j+jb-1)) ? *m : (j+jb-1);
49 /* Apply to columns 1 throughj j-1 */
53 F77_FUNC(dlaswp,DLASWP)(&k,a,lda,&j,&i,ipiv,&l);
55 /* Apply to cols. j+jb through n */
59 F77_FUNC(dlaswp,DLASWP)(&k,&(a[(j+jb-1)*(*lda)+0]),lda,&j,&i,ipiv,&l);
60 /* Compute block row of U */
62 F77_FUNC(dtrsm,DTRSM)("Left","Lower","No transpose","Unit",&jb,&k,&one,
63 &(a[(j-1)*(*lda)+(j-1)]),lda,&(a[(j+jb-1)*(*lda)+(j-1)]),lda);
66 /* Update trailing submatrix */
69 F77_FUNC(dgemm,DGEMM)("No transpose","No transpose",&k,&i,&jb,&minusone,
70 &(a[(j-1)*(*lda)+(j+jb-1)]),lda,
71 &(a[(j+jb-1)*(*lda)+(j-1)]),lda,&one,
72 &(a[(j+jb-1)*(*lda)+(j+jb-1)]),lda);