← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/input/linalg.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/input/linalg.cc
 
1
/*! \file
 
2
    \ingroup INPUT
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdlib>
 
6
#include <cmath>
 
7
#include "linalg.h"
 
8
 
 
9
namespace psi { namespace input {
 
10
 
 
11
FLOAT** create_matrix(int a, int b)
 
12
{
 
13
  int i;
 
14
  FLOAT** M;
 
15
 
 
16
  if (a>=0 && b>=0) {
 
17
    M = (FLOAT**) malloc(sizeof(FLOAT*)*a);
 
18
    M[0] = (FLOAT *) malloc(sizeof(FLOAT)*a*b);
 
19
    for(i=1; i<a; i++)
 
20
      M[i] = M[i-1] + b;
 
21
  }
 
22
  else
 
23
    M = NULL;
 
24
 
 
25
  return M;
 
26
}
 
27
 
 
28
void delete_matrix(FLOAT** M)
 
29
{
 
30
  if (M) {
 
31
    free(M[0]);
 
32
    free(M);
 
33
  }
 
34
}
 
35
 
 
36
void print_matrix(FLOAT** a, int m, int n, FILE* out)
 
37
{
 
38
  int ii,jj,kk,nn,ll;
 
39
  int i,j,k;
 
40
 
 
41
  ii=0;jj=0;
 
42
L200:
 
43
  ii++;
 
44
  jj++;
 
45
  kk=10*jj;
 
46
  nn=n;
 
47
  if (nn > kk) nn=kk;
 
48
  ll = 2*(nn-ii+1)+1;
 
49
  fprintf (out,"\n");
 
50
  for (i=ii; i <= nn; i++) fprintf(out,"       %5d",i);
 
51
  fprintf (out,"\n");
 
52
  for (i=0; i < m; i++) {
 
53
    fprintf (out,"\n%5d",i+1);
 
54
    for (j=ii-1; j < nn; j++) {
 
55
#if LONG_DOUBLE
 
56
      fprintf (out,"%12.7Lf",a[i][j]);
 
57
#else
 
58
      fprintf (out,"%12.7lf",a[i][j]);
 
59
#endif
 
60
    }
 
61
  }
 
62
  fprintf (out,"\n");
 
63
  if (n <= kk) {
 
64
    fflush(out);
 
65
    return;
 
66
  }
 
67
  ii=kk; goto L200;
 
68
}
 
69
 
 
70
FLOAT** convert_matrix(double **m, int a, int b, int transpose)
 
71
{
 
72
  int nrow = transpose ? b : a;
 
73
  int ncol = transpose ? a : b;
 
74
  int nelem, elem, row, col;
 
75
  FLOAT* Melem;
 
76
  double* melem;
 
77
 
 
78
  FLOAT** M = create_matrix(nrow, ncol);
 
79
  if (!transpose) {
 
80
    Melem = M[0];
 
81
    for(row=0; row<nrow; row++)
 
82
      for(col=0; col<ncol; col++) {
 
83
        (*Melem) = (FLOAT) m[row][col];
 
84
        ++Melem;
 
85
      }
 
86
  }
 
87
  else {
 
88
    Melem = M[0];
 
89
    for(row=0; row<nrow; row++)
 
90
      for(col=0; col<ncol; col++) {
 
91
        (*Melem) = (FLOAT) m[col][row];
 
92
        ++Melem;
 
93
      }
 
94
  }
 
95
 
 
96
  return M;
 
97
}
 
98
 
 
99
int matrix_mult(FLOAT** A, int arow, int acol, FLOAT** B, int brow, int bcol, FLOAT** C)
 
100
{
 
101
  int nlink = acol;
 
102
  int row, col, link;
 
103
  FLOAT *aelem, *belem;
 
104
  FLOAT tmp;
 
105
  if (acol != brow)
 
106
    return 1;
 
107
 
 
108
  for(row=0; row<arow; row++) {
 
109
    for(col=0; col<bcol; col++) {
 
110
      aelem = A[row];
 
111
      belem = B[0] + col;
 
112
      tmp = 0.0;
 
113
      for(link=0; link<nlink; link++) {
 
114
        tmp += (*aelem) * (*belem);
 
115
        ++aelem;
 
116
        belem += bcol;
 
117
      }
 
118
      C[row][col] = tmp;
 
119
    }
 
120
  }
 
121
 
 
122
  return 0;
 
123
}
 
124
 
 
125
 
 
126
#define TINY 1.0E-20
 
127
 
 
128
void lu_decom(FLOAT** a, int n, int* indx, FLOAT* d)
 
129
{
 
130
  int i,imax,j,k;
 
131
  FLOAT big,dum,sum,temp;
 
132
  FLOAT* vv = (FLOAT *) malloc(n*sizeof(FLOAT));
 
133
 
 
134
  *d = 1.0;
 
135
 
 
136
  for (i=0; i < n ; i++) {
 
137
    big=0.0;
 
138
    for (j=0; j < n; j++) {
 
139
      if ((temp=FABS(a[i][j])) > big) big=temp;
 
140
    }
 
141
    if (big == 0.0) {
 
142
      *d = 0.0;
 
143
      return;
 
144
    }
 
145
    vv[i] = 1.0/big;
 
146
  }
 
147
  for (j=0; j < n ; j++) {
 
148
    for (i=0; i < j ; i++) {
 
149
      sum = a[i][j];
 
150
      for (k=0; k < i ; k++) sum -= a[i][k]*a[k][j];
 
151
      a[i][j] = sum;
 
152
    }
 
153
    big = 0.0;
 
154
    for (i=j ; i < n ; i++) {
 
155
      sum=a[i][j];
 
156
      for (k=0; k < j ; k++) sum -= a[i][k]*a[k][j];
 
157
      a[i][j] = sum;
 
158
      if ((dum=vv[i]*fabs(sum)) >= big) {
 
159
        big = dum;
 
160
        imax = i;
 
161
      }
 
162
    }
 
163
    if (j != imax) {
 
164
      for (k=0; k < n; k++) {
 
165
        dum=a[imax][k];
 
166
        a[imax][k]=a[j][k];
 
167
        a[j][k]=dum;
 
168
      }
 
169
      *d = -(*d);
 
170
      vv[imax]=vv[j];
 
171
    }
 
172
    indx[j]=imax;
 
173
    if (a[j][j] == 0.0) a[j][j] = TINY;
 
174
    if (j != n-1) {
 
175
      dum = 1.0/a[j][j];
 
176
      for (i=j+1; i < n ; i++) a[i][j] *= dum;
 
177
    }
 
178
  }
 
179
  free(vv);
 
180
}
 
181
 
 
182
}} // namespace psi::input