ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/openbabel/matrix.cpp
Revision: 2440
Committed: Wed Nov 16 19:42:11 2005 UTC (18 years, 8 months ago) by tim
File size: 11416 byte(s)
Log Message:
adding openbabel

File Contents

# User Rev Content
1 tim 2440 /**********************************************************************
2     matrix.cpp - Operations on arbitrary-sized matrix.
3    
4     Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5     Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6    
7     This file is part of the Open Babel project.
8     For more information, see <http://openbabel.sourceforge.net/>
9    
10     This program is free software; you can redistribute it and/or modify
11     it under the terms of the GNU General Public License as published by
12     the Free Software Foundation version 2 of the License.
13    
14     This program 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
17     GNU General Public License for more details.
18     ***********************************************************************/
19    
20     #include "matrix.hpp"
21     #include "vector3.hpp"
22    
23     using namespace std;
24    
25     namespace OpenBabel
26     {
27    
28     void print_matrix(std::vector<std::vector<double> > &m)
29     {
30     unsigned int i,j;
31    
32     for (i = 0; i < m.size(); i++)
33     {
34     for (j = 0; j < m[i].size(); j++)
35     printf("%5.2f",m[i][j]);
36     printf("\n");
37     }
38     }
39    
40     void print_matrix_f(double *m, int rows, int cols)
41     {
42     int i,j,idx;
43    
44     for (i = 0; i < rows; i++)
45     {
46     idx = i * cols;
47     for (j = 0; j < cols; j++)
48     printf("%5.2f",m[idx+j]);
49     printf("\n");
50     }
51     }
52    
53     void print_matrix_ff(double **m, int rows, int cols)
54     {
55     int i,j;
56    
57     for (i = 0; i < rows; i++)
58     {
59     for (j = 0; j < cols; j++)
60     printf("%5.2f",m[i][j]);
61     printf("\n");
62     }
63     }
64    
65     bool mult_matrix(std::vector<std::vector<double> > &c,
66     std::vector<std::vector<double> > &a,
67     std::vector<std::vector<double> > &b)
68     {
69     unsigned int i,j,k;
70    
71     if (a.size() != b.size())
72     return(false);
73    
74     c.resize(a.size());
75    
76     for (i = 0; i < a.size(); i++)
77     {
78     c[i].resize(b[i].size());
79     for (j = 0; j < b[i].size(); j++)
80     {
81     c[i][j] = 0.0;
82     for (k = 0; k < a[i].size(); k++)
83     c[i][j] = c[i][j] + a[i][k] * b[k][j];
84     }
85     }
86    
87     return(true);
88     }
89    
90     bool mult_matrix_f(double *c, double *a, double *b, int rows, int cols)
91     {
92     int i,j,k,idx;
93    
94     for ( i = 0 ; i < rows ; i++ )
95     {
96     idx = i * cols;
97     for ( j = 0; j < cols ; j++ )
98     {
99     c[idx+j] = 0.0;
100     for ( k = 0; k < cols ; k++ )
101     c[idx+j] = c[idx+j] + a[idx+k] * b[(k*cols)+j];
102     }
103     }
104    
105     return(true);
106     }
107    
108     bool mult_matrix_ff(double **c, double **a, double **b, int rows, int cols)
109     {
110     int i,j,k;
111    
112     for ( i = 0 ; i < rows ; i++ )
113     for ( j = 0; j < cols ; j++ )
114     {
115     c[i][j] = 0.0;
116     for ( k = 0; k < cols ; k++ )
117     c[i][j] = c[i][j] + a[i][k] * b[k][j];
118     }
119    
120     return(true);
121     }
122    
123     bool invert_matrix(std::vector<std::vector<double> > &mat, double &det)
124     {
125     int i, j, k, m, n, row = 0, col = 0;
126     double tempo, big, pvt;
127    
128     vector<int> pvt_ind;
129     vector<vector<int> > index;
130    
131     int cols = mat[0].size();
132     int rows = mat.size();
133    
134     pvt_ind.resize(mat[0].size());
135    
136     index.resize(mat.size());
137     for (i = 0; (unsigned)i < mat.size(); i++)
138     index[i].resize(2);
139    
140     // make sure we have a square matrix
141     // #rows == #cols;
142     if (cols != rows)
143     {
144     det = 0.0;
145     return(false);
146     }
147    
148     det = 1.0;
149    
150     for (i = 0; i < cols; i++)
151     pvt_ind[i] = rows+1;
152    
153     for (i = 0; i < cols; i++)
154     {
155     big = 0.0;
156     for (j = 0; j < cols; j++)
157     {
158     if (pvt_ind[j] != 0)
159     for (k = 0; k < cols; k++)
160     {
161     if (fabs(big) < fabs(mat[j][k]))
162     {
163     row = j;
164     col = k;
165     big = mat[j][k];
166     }
167     }
168     }
169    
170     pvt_ind[col]++;
171     if (row != col)
172     {
173     det = -det;
174     for (m = 0; m < cols; m++)
175     {
176     tempo = mat[row][m];
177     mat[row][m] = mat[col][m];
178     mat[col][m] = tempo;
179     }
180     }
181    
182     index[i][0] = row;
183     index[i][1] = col;
184     pvt = mat[col][col];
185     det *= pvt;
186    
187     mat[col][col] = 1.0;
188    
189     for (m = 0; m < cols; m++)
190     mat[col][m] /= pvt;
191    
192     for (n = 0; n < cols; n++)
193     if (n != col)
194     {
195     tempo = mat[n][col];
196     mat[n][col] = 0.0;
197     for (m = 0; m < cols; m++)
198     mat[n][m] -= mat[col][m] * tempo;
199     }
200     }
201    
202     for (i = 0; i < cols; i++)
203     {
204     m = cols - 1;
205     if (index[m][0] != index[m][1])
206     {
207     row = index[m][0];
208     col = index[m][1];
209     for (k = 0; k < cols; k++)
210     {
211     tempo = mat[k][row];
212     mat[k][row] = mat[k][col];
213     mat[k][col] = tempo;
214     }
215     }
216     }
217    
218     return(true);
219     }
220    
221     bool invert_matrix_f(double *mat, double &det, int rows, int cols)
222     {
223     int i, j, k, m, n, row = 0, col = 0, idx1, idx2;
224     double tempo, big, pvt;
225    
226     vector<int> pvt_ind;
227     vector<vector<int> > index;
228    
229     pvt_ind.resize(cols);
230     index.resize(rows);
231    
232     for (i = 0; i < rows; i++)
233     index[i].resize(2);
234    
235     // make sure we have a square matrix
236     // #rows == #cols;
237     if (cols != rows)
238     {
239     det = 0.0;
240     return(false);
241     }
242    
243     det = 1.0;
244    
245     for (i = 0; i < cols; i++)
246     pvt_ind[i] = rows+1;
247    
248     for (i = 0; i < cols; i++)
249     {
250     big = 0.0;
251     for (j = 0; j < cols; j++)
252     {
253     if (pvt_ind[j] != 0)
254     {
255     idx1 = (j * cols);
256     for (k = 0; k < cols; k++)
257     {
258     idx2 = idx1 + k;
259     if (fabs(big) < fabs(mat[idx2]))
260     {
261     row = j;
262     col = k;
263     big = mat[idx2];
264     }
265     }
266     }
267     }
268    
269     pvt_ind[col]++;
270     if (row != col)
271     {
272     det = -det;
273     idx1 = row * cols;
274     idx2 = col * cols;
275     for (m = 0; m < cols; m++)
276     {
277     tempo = mat[idx1+m];
278     mat[idx1+m] = mat[idx2+m];
279     mat[idx2+m] = tempo;
280     }
281     }
282    
283     index[i][0] = row;
284     index[i][1] = col;
285    
286     idx1 = (col*cols);
287     pvt = mat[idx1+col];
288     det *= pvt;
289    
290     mat[idx1+col] = 1.0;
291    
292     for (m = 0; m < cols; m++)
293     mat[idx1+m] /= pvt;
294    
295     for (n = 0; n < cols; n++)
296     if (n != col)
297     {
298     idx1 = n * cols;
299     tempo = mat[idx1 + col];
300     mat[idx1 + col] = 0.0;
301    
302     idx2 = col * cols;
303     for (m = 0; m < cols; m++)
304     mat[idx1 + m] -= mat[idx2 + m] * tempo;
305     }
306     }
307    
308     for (i = 0; i < cols; i++)
309     {
310     m = cols - 1;
311     if (index[m][0] != index[m][1])
312     {
313     row = index[m][0];
314     col = index[m][1];
315     for (k = 0; k < cols; k++)
316     {
317     idx1 = (k * cols);
318     idx2 = idx1 + col;
319     idx1 += row;
320    
321     tempo = mat[idx1];
322     mat[idx1] = mat[idx2];
323     mat[idx2] = tempo;
324     }
325     }
326     }
327    
328     return(true);
329     }
330    
331     bool invert_matrix_ff(double **mat, double &det, int rows, int cols)
332     {
333     int i, j, k, m, n, row = 0, col = 0;
334     double tempo, big, pvt;
335    
336     vector<int> pvt_ind;
337     vector<vector<int> > index;
338    
339     pvt_ind.resize(cols);
340     index.resize(rows);
341    
342     for (i = 0; i < rows; i++)
343     index[i].resize(2);
344    
345     // make sure we have a square matrix
346     // #rows == #cols;
347     if (cols != rows)
348     {
349     det = 0.0;
350     return(false);
351     }
352    
353     det = 1.0;
354    
355     for (i = 0; i < cols; i++)
356     pvt_ind[i] = rows+1;
357    
358     for (i = 0; i < cols; i++)
359     {
360     big = 0.0;
361     for (j = 0; j < cols; j++)
362     {
363     if (pvt_ind[j] != 0)
364     for (k = 0; k < cols; k++)
365     {
366     if (fabs(big) < fabs(mat[j][k]))
367     {
368     row = j;
369     col = k;
370     big = mat[j][k];
371     }
372     }
373     }
374    
375     pvt_ind[col]++;
376     if (row != col)
377     {
378     det = -det;
379     for (m = 0; m < cols; m++)
380     {
381     tempo = mat[row][m];
382     mat[row][m] = mat[col][m];
383     mat[col][m] = tempo;
384     }
385     }
386    
387     index[i][0] = row;
388     index[i][1] = col;
389     pvt = mat[col][col];
390     det *= pvt;
391    
392     mat[col][col] = 1.0;
393    
394     for (m = 0; m < cols; m++)
395     mat[col][m] /= pvt;
396    
397     for (n = 0; n < cols; n++)
398     if (n != col)
399     {
400     tempo = mat[n][col];
401     mat[n][col] = 0.0;
402     for (m = 0; m < cols; m++)
403     mat[n][m] -= mat[col][m] * tempo;
404     }
405     }
406    
407     for (i = 0; i < cols; i++)
408     {
409     m = cols - 1;
410     if (index[m][0] != index[m][1])
411     {
412     row = index[m][0];
413     col = index[m][1];
414     for (k = 0; k < cols; k++)
415     {
416     tempo = mat[k][row];
417     mat[k][row] = mat[k][col];
418     mat[k][col] = tempo;
419     }
420     }
421     }
422    
423     return(true);
424     }
425    
426     bool convert_matrix_f(std::vector<std::vector<double> > &src, double *dst)
427     {
428     unsigned int i, j, idx = 0;
429    
430     for ( i = 0 ; i < src.size() ; i++ )
431     for ( j = 0 ; j < src[i].size() ; j++ )
432     dst[idx++] = src[i][j];
433    
434     return true;
435     }
436    
437     bool convert_matrix_ff(std::vector<std::vector<double> > &src, double **dst)
438     {
439     unsigned int i, j;
440    
441     for ( i = 0 ; i < src.size() ; i++ )
442     for ( j = 0 ; j < src[i].size() ; j++ )
443     dst[i][j] = src[i][j];
444    
445     return true;
446     }
447    
448     bool convert_matrix_f(double *src, std::vector<std::vector<double> > &dst,
449     int rows, int cols)
450     {
451     int i, j, idx;
452    
453     dst.resize(rows);
454     for ( i = 0 ; i < rows ; i++ )
455     {
456     idx = i * cols;
457     dst[i].resize(cols);
458     for ( j = 0 ; j < cols ; j++ )
459     dst[i][j] = src[idx+j];
460     }
461    
462     return true;
463     }
464    
465     bool convert_matrix_ff(double **src, std::vector<std::vector<double> > &dst,
466     int rows, int cols)
467     {
468     int i, j;
469    
470     dst.resize(rows);
471     for ( i = 0 ; i < rows ; i++ )
472     {
473     dst[i].resize(cols);
474     for ( j = 0 ; j < cols ; j++ )
475     dst[i][j] = src[i][j];
476     }
477    
478     return true;
479     }
480    
481     bool convert_matrix_f_ff(double *src, double **dst, int rows, int cols)
482     {
483     int i, j, idx;
484    
485     for ( i = 0 ; i < rows ; i++ )
486     {
487     idx = i * cols;
488     for ( j = 0 ; j < cols ; j++ )
489     dst[i][j] = src[idx+j];
490     }
491    
492     return true;
493     }
494    
495     bool convert_matrix_ff_f(double **src, double *dst, int rows, int cols)
496     {
497     int i, j, idx;
498    
499     for ( i = 0 ; i < rows ; i++ )
500     {
501     idx = i * cols;
502     for ( j = 0 ; j < cols ; j++ )
503     dst[idx+j] = src[i][j];
504     }
505    
506     return true;
507     }
508    
509     } // end namespace OpenBabel
510    
511     //! \file matrix.cpp
512     //! \brief Operations on arbitrary-sized matrix.
513