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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathMatv_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/mathMatv_invert.c,v $
 37  * Date    :  $Date: 2003/09/23 11:32:20 $
 38  * Version :  $Revision: 1.4 $
 39  * CVS Id  :  $Id: mathMatv_invert.c,v 1.4 2003/09/23 11:32:20 matts Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Matrix inversion.
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief   Matrix inversion.  
 50  *
 51  *  NB, the mat_invert() function takes a Mat structure as an argument, NOT a Matrix.
 52  *  The invert() function simply inverts an array of doubles.
 53 */
 54 
 55 #include "mathMatv_invert.h"
 56 
 57 #if HAVE_CONFIG_H
 58 #include <config.h>
 59 #endif
 60 
 61 #include <math.h>
 62 #include <tina/sys/sysDef.h>
 63 #include <tina/sys/sysPro.h>
 64 #include <tina/math/math_MatrDef.h>
 65 
 66 static void     invert(int n, double **a)
 67 {
 68         int             i, j, k;
 69         int            *p = tvector_alloc(0, n, int);
 70         double          h, q, s, sup, pivot;
 71         for (k = 0; k < n; k++)
 72         {
 73                 sup = 0.0;
 74                 p[k] = 0;
 75                 for (i = k; i < n; i++)
 76                 {
 77                         s = 0.0;
 78                         for (j = k; j < n; j++)
 79                                 s += fabs(a[i][j]);
 80                         q = fabs(a[i][k]) / s;
 81                         if (sup < q)
 82                         {
 83                                 sup = q;
 84                                 p[k] = i;
 85                         }
 86                 }
 87                 if (sup == 0.0)
 88                         return;
 89                 if (p[k] != k)
 90                         for (j = 0; j < n; j++)
 91                         {
 92                                 h = a[k][j];
 93                                 a[k][j] = a[p[k]][j];
 94                                 a[p[k]][j] = h;
 95                         }
 96                 pivot = a[k][k];
 97                 for (j = 0; j < n; j++)
 98                         if (j != k)
 99                         {
100                                 a[k][j] = -a[k][j] / pivot;
101                                 for (i = 0; i < n; i++)
102                                         if (i != k)
103                                                 a[i][j] += a[i][k] * a[k][j];
104                         }
105                 for (i = 0; i < n; i++)
106                         a[i][k] = a[i][k] / pivot;
107                 a[k][k] = 1.0 / pivot;
108         }
109         for (k = n - 1; k >= 0; k--)
110                 if (p[k] != k)
111                         for (i = 0; i < n; i++)
112                         {
113                                 h = a[i][k];
114                                 a[i][k] = a[i][p[k]];
115                                 a[i][p[k]] = h;
116                         }
117 
118         tvector_free(p, 0, int);
119 }
120 
121 /*
122 Invert matrix a in place using pivoting method above.
123 */
124 void            mat_invert(Mat * a)
125 {
126         if (a == NULL)
127                 return;
128         invert(a->n, a->el);
129 }
130 

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