Fixing copyright issues and code contributors
[alexxy/gromacs.git] / src / gmxlib / gmx_lapack / sgetrf.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2012,2013, by the GROMACS development team, led by
5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6  * others, as listed in the AUTHORS file in the top-level source
7  * directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 #include "gmx_blas.h"
36 #include "gmx_lapack.h"
37 #include "lapack_limits.h"
38
39 void
40 F77_FUNC(sgetrf,SGETRF)(int *m,
41         int *n,
42         float *a,
43         int *lda,
44         int *ipiv,
45         int *info)
46 {
47   int mindim,jb;
48   int i,j,k,l;
49   int iinfo;
50   float minusone = -1.0;
51   float one = 1.0;
52
53   if(*m<=0 || *n<=0)
54     return;
55
56   *info = 0;
57
58   mindim = (*m < *n) ? *m : *n;
59
60   if(DGETRF_BLOCKSIZE>=mindim) {
61
62     /* unblocked code */
63     F77_FUNC(sgetf2,SGETF2)(m,n,a,lda,ipiv,info);
64
65   } else {
66
67     /* blocked case */
68
69     for(j=1;j<=mindim;j+=DGETRF_BLOCKSIZE) {
70       jb = ( DGETRF_BLOCKSIZE < (mindim-j+1)) ? DGETRF_BLOCKSIZE : (mindim-j+1);
71       /* factor diag. and subdiag blocks and test for singularity */
72       k = *m-j+1;
73       F77_FUNC(sgetf2,SGETF2)(&k,&jb,&(a[(j-1)*(*lda)+(j-1)]),lda,&(ipiv[j-1]),&iinfo);
74       
75       if(*info==0 && iinfo>0)
76         *info = iinfo + j - 1;
77
78       /* adjust pivot indices */
79       k = (*m < (j+jb-1)) ? *m : (j+jb-1);
80       for(i=j;i<=k;i++)
81         ipiv[i-1] += j - 1;
82
83       /* Apply to columns 1 throughj j-1 */
84       k = j - 1;
85       i = j + jb - 1;
86       l = 1;
87       F77_FUNC(slaswp,SLASWP)(&k,a,lda,&j,&i,ipiv,&l);
88       if((j+jb)<=*n) {
89         /* Apply to cols. j+jb through n */
90         k = *n-j-jb+1;
91         i = j+jb-1;
92         l = 1;
93         F77_FUNC(slaswp,SLASWP)(&k,&(a[(j+jb-1)*(*lda)+0]),lda,&j,&i,ipiv,&l);
94         /* Compute block row of U */
95         k = *n-j-jb+1;
96         F77_FUNC(strsm,STRSM)("Left","Lower","No transpose","Unit",&jb,&k,&one,
97                &(a[(j-1)*(*lda)+(j-1)]),lda,&(a[(j+jb-1)*(*lda)+(j-1)]),lda);
98
99         if((j+jb)<=*m) {
100           /* Update trailing submatrix */
101           k = *m-j-jb+1;
102           i = *n-j-jb+1;
103           F77_FUNC(sgemm,SGEMM)("No transpose","No transpose",&k,&i,&jb,&minusone,
104                  &(a[(j-1)*(*lda)+(j+jb-1)]),lda,
105                  &(a[(j+jb-1)*(*lda)+(j-1)]),lda,&one,
106                  &(a[(j+jb-1)*(*lda)+(j+jb-1)]),lda);
107         }
108
109       }
110     }
111   }
112 }