Back to index

radiance  4R0+20100331
invmat4.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: invmat4.c,v 2.3 2003/02/25 02:47:21 greg Exp $";
00003 #endif
00004 /*
00005  * invmat4 - computes the inverse of mat into inverse.  Returns 1
00006  * if there exists an inverse, 0 otherwise.  It uses Gaussian Elimination
00007  * method.
00008  */
00009 
00010 #include "copyright.h"
00011 
00012 #include "mat4.h"
00013 
00014 
00015 #define  ABS(x)             ((x)>=0 ? (x) : -(x))
00016 
00017 #define SWAP(a,b,t) (t=a,a=b,b=t)
00018 
00019 
00020 int
00021 invmat4(inverse, mat)
00022 MAT4  inverse, mat;
00023 {
00024        MAT4  m4tmp;
00025        register int i,j,k;
00026        register double temp;
00027 
00028        copymat4(m4tmp, mat);
00029                                    /* set inverse to identity */
00030        setident4(inverse);
00031 
00032        for(i = 0; i < 4; i++) {
00033               /* Look for row with largest pivot and swap rows */
00034               temp = FTINY; j = -1;
00035               for(k = i; k < 4; k++)
00036                      if(ABS(m4tmp[k][i]) > temp) {
00037                             temp = ABS(m4tmp[k][i]);
00038                             j = k;
00039                             }
00040               if(j == -1)   /* No replacing row -> no inverse */
00041                      return(0);
00042               if (j != i)
00043                      for(k = 0; k < 4; k++) {
00044                             SWAP(m4tmp[i][k],m4tmp[j][k],temp);
00045                             SWAP(inverse[i][k],inverse[j][k],temp);
00046                             }
00047 
00048               temp = m4tmp[i][i];
00049               for(k = 0; k < 4; k++) {
00050                      m4tmp[i][k] /= temp;
00051                      inverse[i][k] /= temp;
00052                      }
00053               for(j = 0; j < 4; j++) {
00054                      if(j != i) {
00055                             temp = m4tmp[j][i];
00056                             for(k = 0; k < 4; k++) {
00057                                    m4tmp[j][k] -= m4tmp[i][k]*temp;
00058                                    inverse[j][k] -= inverse[i][k]*temp;
00059                                    }
00060                             }
00061                      }
00062               }
00063        return(1);
00064 }