~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/matrix/mat_invert.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)Matrix invertion (various types & shapes)
  2   @(#)(arithmetic always done in double precision)
  3  */
  4 
  5 #include <math.h>
  6 #include <tina/sys.h>
  7 #include <tina/sysfuncs.h>
  8 #include <tina/math.h>
  9 #include <tina/mathfuncs.h>
 10 
 11 int     matrix_get();
 12 double  matrix_getf();
 13 
 14 void    matrix_put();
 15 void    matrix_putf();
 16 
 17 double **darray_invert(double **a1, int n);
 18 double **dupper_invert(double **u, int n);
 19 double **dlower_invert(double **l, int n);
 20 
 21 Matrix *matrix_alloc();
 22 Matrix *matrix_build();
 23 Matrix *matrix_copy();
 24 Matrix *matrix_cast();
 25 Matrix *matrix_fill();
 26 Matrix *dmatrix_invert(Matrix * mat);
 27 
 28 Vartype matrix_sup_vtype();
 29 
 30 Matrix *matrix_invert(Matrix * mat)
 31 {
 32     Matrix *inverse;
 33     Matrix *inverse1;
 34     Matrix *m;
 35 
 36     if (mat == NULL)
 37         return (mat);
 38 
 39     /** make copy of type double if necessary **/
 40 
 41     if (mat->vtype != double_v)
 42         m = matrix_cast(mat, double_v);
 43     else
 44         m = mat;
 45 
 46     inverse1 = dmatrix_invert(m);
 47 
 48     switch (mat->vtype)
 49     {
 50     case int_v:
 51     case float_v:
 52         inverse = matrix_cast(inverse1, float_v);
 53         matrix_free((Matrix *) inverse1);
 54         break;
 55     case double_v:
 56         inverse = inverse1;
 57         break;
 58     }
 59 
 60     /** free copy if necessary **/
 61 
 62     if (mat->vtype != double_v)
 63         matrix_free((Matrix *) m);
 64 
 65     return (inverse);
 66 }
 67 
 68 Matrix *dmatrix_invert(Matrix * mat)
 69 {
 70     Matrix *inverse;
 71     double **el;
 72     int     m;
 73 
 74     m = MIN(mat->m, mat->n);
 75     switch (mat->shape)
 76     {
 77     case matrix_full:
 78         el = darray_invert(mat->el.double_v, m);
 79         inverse = matrix_build(m, m, matrix_full, double_v, (void *) el);
 80         break;
 81     case matrix_lower:
 82         el = dlower_invert(mat->el.double_v, m);
 83         inverse = matrix_build(m, m, matrix_upper, double_v, (void *) el);
 84         break;
 85     case matrix_upper:
 86         el = dupper_invert(mat->el.double_v, m);
 87         inverse = matrix_build(m, m, matrix_lower, double_v, (void *) el);
 88         break;
 89     default:
 90         {
 91             Matrix *mat1 = matrix_fill(mat);
 92 
 93             inverse = dmatrix_invert(mat1);
 94             matrix_free((Matrix *) mat1);
 95             break;
 96         }
 97     }
 98     return (inverse);
 99 }
100 
101 double **darray_invert(double **a1, int n)
102 {
103     int     i, j, k;
104     int    *p = ivector_alloc(0, n);
105     double **a = darray_alloc(0, 0, n, n);
106     double  h, q, s, sup, pivot;
107 
108 
109     for (i = 0; i < n; i++)
110         for (j = 0; j < n; j++)
111             a[i][j] = a1[i][j];
112 
113     for (k = 0; k < n; k++)
114     {
115         sup = 0.0;
116         p[k] = 0;
117         for (i = k; i < n; i++)
118         {
119             s = 0.0;
120             for (j = k; j < n; j++)
121                 s += fabs(a[i][j]);
122             q = fabs(a[i][k]) / s;
123             if (sup < q)
124             {
125                 sup = q;
126                 p[k] = i;
127             }
128         }
129         if (sup == 0.0)
130             return (NULL);
131         if (p[k] != k)
132             for (j = 0; j < n; j++)
133             {
134                 h = a[k][j];
135                 a[k][j] = a[p[k]][j];
136                 a[p[k]][j] = h;
137             }
138         pivot = a[k][k];
139         for (j = 0; j < n; j++)
140             if (j != k)
141             {
142                 a[k][j] = -a[k][j] / pivot;
143                 for (i = 0; i < n; i++)
144                     if (i != k)
145                         a[i][j] += a[i][k] * a[k][j];
146             }
147         for (i = 0; i < n; i++)
148             a[i][k] = a[i][k] / pivot;
149         a[k][k] = 1.0 / pivot;
150     }
151     for (k = n - 1; k >= 0; k--)
152         if (p[k] != k)
153             for (i = 0; i < n; i++)
154             {
155                 h = a[i][k];
156                 a[i][k] = a[i][p[k]];
157                 a[i][p[k]] = h;
158             }
159 
160     ivector_free((void *) p, 0);
161     return (a);
162 }
163 
164 double **dlower_invert(double **l, int n)
165 {
166     int     i, j, k;
167     double  sum;
168     double **u = dupper_alloc(0, n);
169 
170     for (j = n - 1; j >= 0; --j)
171     {
172         u[j][j] = -1.0 / l[j][j];
173         for (i = j + 1; i < n; ++i)
174         {
175             sum = 0.0;
176             for (k = j; k < i; ++k)
177                 sum += l[i][k] * u[k][j];
178             u[i][j] = u[j][j] * sum;
179         }
180     }
181 
182     return (u);
183 }
184 
185 double **dupper_invert(double **u, int n)
186 {
187     int     i, j, k;
188     double  sum;
189     double **l = dlower_alloc(0, n);
190 
191     for (j = 0; j < n; ++j)
192     {
193         l[j][j] = -1.0 / u[j][j];
194         for (i = j - 1; i >= 0; --i)
195         {
196             sum = 0.0;
197             for (k = i; k < j; ++k)
198                 sum += u[i][k] * l[k][j];
199             l[i][j] = l[j][j] * sum;
200         }
201     }
202     return (l);
203 }
204 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.