1 /* 2 3 -Procedure invert_c ( Invert a 3x3 matrix ) 4 5 -Abstract 6 7 Generate the inverse of a 3x3 matrix. 8 9 -Disclaimer 10 11 THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE 12 CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. 13 GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE 14 ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE 15 PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" 16 TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY 17 WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A 18 PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC 19 SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE 20 SOFTWARE AND RELATED MATERIALS, HOWEVER USED. 21 22 IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA 23 BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT 24 LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, 25 INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, 26 REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE 27 REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. 28 29 RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF 30 THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY 31 CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE 32 ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. 33 34 -Required_Reading 35 36 None. 37 38 -Keywords 39 40 MATRIX, MATH 41 42 */ 43 #include <math.h> 44 #include "SpiceUsr.h" 45 #undef invert_c 46 47 invert_c(ConstSpiceDouble m1[3][3],SpiceDouble mout[3][3])48 void invert_c ( ConstSpiceDouble m1 [3][3], 49 SpiceDouble mout[3][3] ) 50 51 /* 52 53 -Brief_I/O 54 55 VARIABLE I/O DESCRIPTION 56 -------- --- -------------------------------------------------- 57 m1 I Matrix to be inverted. 58 mout O Inverted matrix (m1)**-1. If m1 is singular, then 59 mout will be the zero matrix. mout can 60 overwrite m1. 61 62 -Detailed_Input 63 64 m1 An arbitrary 3x3 matrix. The limits on the size of 65 elements of m1 are determined by the process of calculating 66 the cofactors of each element of the matrix. For a 3x3 67 matrix this amounts to the differencing of two terms, each 68 of which consists of the multiplication of two matrix 69 elements. This multiplication must not exceed the range 70 of double precision numbers or else an overflow error will 71 occur. 72 73 -Detailed_Output 74 75 mout is the inverse of m1 and is calculated explicitly using 76 the matrix of cofactors. mout is set to be the zero matrix 77 if m1 is singular. 78 79 -Parameters 80 81 None. 82 83 -Exceptions 84 85 1) No internal checking on the input matrix m1 is performed except on 86 the size of its determinant. Thus it is possible to generate a 87 floating point overflow or underflow in the process of 88 calculating the matrix of cofactors. 89 90 2) If the determinant is less than 10**-16, the matrix is deemed to 91 be singular and the output matrix is filled with zeros. 92 93 -Particulars 94 95 A temporary matrix is used to compute the result, so the output 96 matrix may overwrite the input matrix. 97 98 -Examples 99 100 Suppose that m1 is given by the following matrix equation: 101 102 | 0 -1 0 | 103 m1 = | 0.5 0 0 | 104 | 0 0 1 | 105 106 If invert_c is called as shown 107 108 invert_c (m1, m1); 109 110 then m1 will be set to be: 111 112 | 0 2 0 | 113 m1 = |-1 0 0 | 114 | 0 0 1 | 115 116 -Restrictions 117 118 The input matrix must be such that generating the cofactors will 119 not cause a floating point overflow or underflow. The 120 strictness of this condition depends, of course, on the computer 121 installation and the resultant maximum and minimum values of 122 double precision numbers. 123 124 -Files 125 126 None 127 128 -Author_and_Institution 129 130 W.M. Owen (JPL) 131 132 -Literature_References 133 134 None 135 136 -Version 137 138 -CSPICE Version 1.0.0, 13-SEP-1999 (NJB) (WMO) 139 140 -Index_Entries 141 142 invert a 3x3_matrix 143 144 -& 145 */ 146 147 { /* Begin invert_c */ 148 149 /* 150 Local constants 151 */ 152 153 #define SINGULAR_DET 1.e-16 154 155 156 157 /* 158 Local variables 159 */ 160 SpiceInt i; 161 162 SpiceDouble invdet; 163 SpiceDouble mdet; 164 SpiceDouble mtemp[3][3]; 165 166 167 /* 168 Find the determinant of m1 and check for singularity. 169 */ 170 171 mdet = det_c(m1); 172 173 if ( fabs(mdet) < SINGULAR_DET ) 174 { 175 176 /* 177 The matrix is considered to be singular. 178 */ 179 180 for ( i = 0; i < 9; i++ ) 181 { 182 *( (SpiceDouble*)mout+i ) = 0.; 183 } 184 185 return; 186 } 187 188 189 /* 190 Get the cofactors of each element of m1. 191 */ 192 mtemp[0][0] = ( m1[1][1]*m1[2][2] - m1[2][1]*m1[1][2] ); 193 mtemp[0][1] = -( m1[0][1]*m1[2][2] - m1[2][1]*m1[0][2] ); 194 mtemp[0][2] = ( m1[0][1]*m1[1][2] - m1[1][1]*m1[0][2] ); 195 mtemp[1][0] = -( m1[1][0]*m1[2][2] - m1[2][0]*m1[1][2] ); 196 mtemp[1][1] = ( m1[0][0]*m1[2][2] - m1[2][0]*m1[0][2] ); 197 mtemp[1][2] = -( m1[0][0]*m1[1][2] - m1[1][0]*m1[0][2] ); 198 mtemp[2][0] = ( m1[1][0]*m1[2][1] - m1[2][0]*m1[1][1] ); 199 mtemp[2][1] = -( m1[0][0]*m1[2][1] - m1[2][0]*m1[0][1] ); 200 mtemp[2][2] = ( m1[0][0]*m1[1][1] - m1[1][0]*m1[0][1] ); 201 202 /* 203 Multiply the cofactor matrix by 1/mdet to obtain the inverse matrix. 204 */ 205 206 invdet = 1. / mdet; 207 208 vsclg_c ( invdet, (SpiceDouble *)mtemp, 9, (SpiceDouble *)mout ); 209 210 211 } /* End invert_c */ 212 213