1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2008, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
43 #include "gmx_fatal.h"
44 #include "gmx_matrix.h"
45 #include "gmx_lapack.h"
47 double **alloc_matrix(int n
,int m
)
52 /* There's always time for more pointer arithmetic! */
53 /* This is necessary in order to be able to work with LAPACK */
61 void free_matrix(double **a
,int n
)
70 void matrix_multiply(FILE *fp
,int n
,int m
,double **x
,double **y
,double **z
)
76 fprintf(fp
,"Multiplying %d x %d matrix with a %d x %d matrix\n",
82 fprintf(fp
," %7g",x
[i
][j
]);
92 z
[i
][j
] += x
[k
][i
]*y
[j
][k
];
97 static void dump_matrix(FILE *fp
,const char *title
,int n
,double **a
)
102 fprintf(fp
,"%s\n",title
);
108 fprintf(fp
," %8.2f",a
[i
][j
]);
112 fprintf(fp
,"Prod a[i][i] = %g\n",d
);
115 int matrix_invert(FILE *fp
,int n
,double **a
)
117 int i
,j
,m
,lda
,*ipiv
,lwork
,info
;
118 double **test
=NULL
,**id
,*work
;
123 fprintf(fp
,"Inverting %d square matrix\n",n
);
124 test
= alloc_matrix(n
,n
);
129 test
[i
][j
] = a
[i
][j
];
132 dump_matrix(fp
,"before inversion",n
,a
);
140 F77_FUNC(dgetrf
,DGETRF
)(&n
,&m
,a
[0],&lda
,ipiv
,&info
);
143 dump_matrix(fp
,"after dgetrf",n
,a
);
147 F77_FUNC(dgetri
,DGETRI
)(&n
,a
[0],&lda
,ipiv
,work
,&lwork
,&info
);
150 dump_matrix(fp
,"after dgetri",n
,a
);
158 id
= alloc_matrix(n
,n
);
159 matrix_multiply(fp
,n
,n
,test
,a
,id
);
160 dump_matrix(fp
,"And here is the product of A and Ainv",n
,id
);
171 double multi_regression(FILE *fp
,int nrow
,double *y
,int ncol
,
172 double **xx
,double *a0
)
175 double ax
,chi2
,**a
,**at
,**ata
,*atx
;
177 a
= alloc_matrix(nrow
,ncol
);
178 at
= alloc_matrix(ncol
,nrow
);
179 ata
= alloc_matrix(ncol
,ncol
);
180 for(i
=0; (i
<nrow
); i
++)
181 for(j
=0; (j
<ncol
); j
++)
182 at
[j
][i
] = a
[i
][j
] = xx
[j
][i
];
183 matrix_multiply(fp
,nrow
,ncol
,a
,at
,ata
);
184 if ((row
= matrix_invert(fp
,ncol
,ata
)) != 0) {
185 gmx_fatal(FARGS
,"Matrix inversion failed. Incorrect row = %d.\nThis probably indicates that you do not have sufficient data points, or that some parameters are linearly dependent.",
190 for(i
=0; (i
<ncol
); i
++)
193 for(j
=0; (j
<nrow
); j
++)
194 atx
[i
] += at
[i
][j
]*y
[j
];
196 for(i
=0; (i
<ncol
); i
++)
199 for(j
=0; (j
<ncol
); j
++)
200 a0
[i
] += ata
[i
][j
]*atx
[j
];
203 for(j
=0; (j
<nrow
); j
++)
206 for(i
=0; (i
<ncol
); i
++)
208 chi2
+= sqr(y
[j
]-ax
);
213 free_matrix(at
,ncol
);
214 free_matrix(ata
,ncol
);