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

Linux Cross Reference
Tina5/tina-libs/tina/math/mathMatr_invert.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathMatr_invert.c,v $
 37  * Date    :  $Date: 2005/01/23 19:10:21 $
 38  * Version :  $Revision: 1.6 $
 39  * CVS Id  :  $Id: mathMatr_invert.c,v 1.6 2005/01/23 19:10:21 paul Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Matrix inversion (various types & shapes)
 44  *         (arithmetic always done in double precision)
 45  *
 46  *********
 47 */
 48 /** 
 49  *  @file
 50  *  @brief  Matrix inversion (various types & shapes).  
 51  *
 52  *  arithmetic always done in double precision
 53  * 
 54 */
 55 
 56 #include "mathMatr_invert.h"
 57 
 58 #if HAVE_CONFIG_H
 59 #include <config.h>
 60 #endif
 61 
 62 #include <math.h>
 63 #include <tina/sys/sysDef.h>
 64 #include <tina/sys/sysPro.h>
 65 #include <tina/math/math_MatrDef.h>
 66 #include <tina/math/mathMatr_get.h>
 67 #include <tina/math/mathMatr_put.h>
 68 #include <tina/math/mathMatr_alloc.h>
 69 #include <tina/math/mathMatr_build.h>
 70 #include <tina/math/mathMatr_copy.h>
 71 #include <tina/math/mathMatr_cast.h>
 72 #include <tina/math/mathMatr_fill.h>
 73 #include <tina/math/mathMatr_util.h>
 74 #include <tina/math/mathMatr_free.h>
 75 
 76 Matrix         *matrix_invert(Matrix * mat)
 77 {
 78         Matrix         *inverse=NULL;
 79         Matrix         *inverse1=NULL;
 80         Matrix         *m=NULL;
 81 
 82         if (mat == NULL)
 83                 return (mat);
 84 
 85         /** make copy of type double if necessary **/
 86 
 87         if (mat->vtype != double_v)
 88                 m = matrix_cast(mat, double_v);
 89         else
 90                 m = mat;
 91 
 92         inverse1 = dmatrix_invert(m);
 93 
 94         switch (mat->vtype)
 95         {
 96         case int_v:
 97         case float_v:
 98                 inverse = matrix_cast(inverse1, float_v);
 99                 matrix_free((Matrix *) inverse1);
100                 break;
101         case double_v:
102                 inverse = inverse1;
103                 break;
104         default:
105                 error("matrix_invert: unsupported type", non_fatal);
106                 break;
107         }
108 
109         /** free copy if necessary **/
110 
111         if (mat->vtype != double_v)
112                 matrix_free((Matrix *) m);
113 
114         return (inverse);
115 }
116 
117 Matrix         *dmatrix_invert(Matrix * mat)
118 {
119         Matrix         *inverse;
120         double        **el;
121         int             m;
122 
123         m = MIN(mat->m, mat->n);
124         switch (mat->shape)
125         {
126         case matrix_full:
127                 el = darray_invert(mat->el.double_v, m);
128                 inverse = matrix_build(m, m, matrix_full, double_v, (void *) el);
129                 break;
130         case matrix_lower:
131                 el = dlower_invert(mat->el.double_v, m);
132                 inverse = matrix_build(m, m, matrix_upper, double_v, (void *) el);
133                 break;
134         case matrix_upper:
135                 el = dupper_invert(mat->el.double_v, m);
136                 inverse = matrix_build(m, m, matrix_lower, double_v, (void *) el);
137                 break;
138         default:
139                 {
140                         Matrix         *mat1 = matrix_fill(mat);
141 
142                         inverse = dmatrix_invert(mat1);
143                         matrix_free((Matrix *) mat1);
144                         break;
145                 }
146         }
147         return (inverse);
148 }
149 
150 double        **darray_invert(double **a1, int n)
151 {
152         int             i, j, k;
153         int            *p = ivector_alloc(0, n);
154         double        **a = darray_alloc(0, 0, n, n);
155         double          h, q, s, sup, pivot;
156 
157 
158         for (i = 0; i < n; i++)
159                 for (j = 0; j < n; j++)
160                         a[i][j] = a1[i][j];
161 
162         for (k = 0; k < n; k++)
163         {
164                 sup = 0.0;
165                 p[k] = 0;
166                 for (i = k; i < n; i++)
167                 {
168                         s = 0.0;
169                         for (j = k; j < n; j++)
170                                 s += fabs(a[i][j]);
171                         q = fabs(a[i][k]) / s;
172                         if (sup < q)
173                         {
174                                 sup = q;
175                                 p[k] = i;
176                         }
177                 }
178                 if (sup == 0.0)
179                         return (NULL);
180                 if (p[k] != k)
181                         for (j = 0; j < n; j++)
182                         {
183                                 h = a[k][j];
184                                 a[k][j] = a[p[k]][j];
185                                 a[p[k]][j] = h;
186                         }
187                 pivot = a[k][k];
188                 for (j = 0; j < n; j++)
189                         if (j != k)
190                         {
191                                 a[k][j] = -a[k][j] / pivot;
192                                 for (i = 0; i < n; i++)
193                                         if (i != k)
194                                                 a[i][j] += a[i][k] * a[k][j];
195                         }
196                 for (i = 0; i < n; i++)
197                         a[i][k] = a[i][k] / pivot;
198                 a[k][k] = 1.0 / pivot;
199         }
200         for (k = n - 1; k >= 0; k--)
201                 if (p[k] != k)
202                         for (i = 0; i < n; i++)
203                         {
204                                 h = a[i][k];
205                                 a[i][k] = a[i][p[k]];
206                                 a[i][p[k]] = h;
207                         }
208 
209         ivector_free(p, 0);
210         return (a);
211 }
212 
213 double        **dlower_invert(double **l, int n)
214 {
215         int             i, j, k;
216         double          sum;
217         double        **u = dupper_alloc(0, n);
218 
219         for (j = n - 1; j >= 0; --j)
220         {
221                 u[j][j] = -1.0 / l[j][j];
222                 for (i = j + 1; i < n; ++i)
223                 {
224                         sum = 0.0;
225                         for (k = j; k < i; ++k)
226                                 sum += l[i][k] * u[k][j];
227                         u[i][j] = u[j][j] * sum;
228                 }
229         }
230 
231         return (u);
232 }
233 
234 double        **dupper_invert(double **u, int n)
235 {
236         int             i, j, k;
237         double          sum;
238         double        **l = dlower_alloc(0, n);
239 
240         for (j = 0; j < n; ++j)
241         {
242                 l[j][j] = -1.0 / u[j][j];
243                 for (i = j - 1; i >= 0; --i)
244                 {
245                         sum = 0.0;
246                         for (k = i; k < j; ++k)
247                                 sum += u[i][k] * l[k][j];
248                         l[i][j] = l[j][j] * sum;
249                 }
250         }
251         return (l);
252 }
253 

~ [ 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.