1 /* 2 3 -Procedure mxmtg_c ( Matrix times matrix transpose, general dimension ) 4 5 -Abstract 6 7 Multiply a matrix and the transpose of a matrix, both of 8 arbitrary size. 9 10 -Disclaimer 11 12 THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE 13 CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. 14 GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE 15 ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE 16 PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" 17 TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY 18 WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A 19 PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC 20 SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE 21 SOFTWARE AND RELATED MATERIALS, HOWEVER USED. 22 23 IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA 24 BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT 25 LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, 26 INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, 27 REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE 28 REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. 29 30 RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF 31 THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY 32 CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE 33 ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. 34 35 -Required_Reading 36 37 None. 38 39 -Keywords 40 41 MATRIX 42 43 */ 44 #include "SpiceUsr.h" 45 #include "SpiceZfc.h" 46 #include "SpiceZst.h" 47 #include "SpiceZmc.h" 48 #include "SpiceZim.h" 49 #undef mxmtg_c 50 51 mxmtg_c(const void * m1,const void * m2,SpiceInt nrow1,SpiceInt nc1c2,SpiceInt nrow2,void * mout)52 void mxmtg_c ( const void * m1, 53 const void * m2, 54 SpiceInt nrow1, 55 SpiceInt nc1c2, 56 SpiceInt nrow2, 57 void * mout ) 58 /* 59 60 -Brief_I/O 61 62 VARIABLE I/O DESCRIPTION 63 -------- --- -------------------------------------------------- 64 m1 I Left-hand matrix to be multiplied. 65 m2 I Right-hand matrix whose transpose is to be multiplied 66 nrow1 I Row dimension of m1 and row dimension of mout. 67 nc1c2 I Column dimension of m1 and column dimension of m2. 68 nrow2 I Row dimension of m2 and column dimension of mout. 69 mout O Product matrix. 70 71 -Detailed_Input 72 73 m1 may be any double precision matrix of arbitrary size. 74 75 m2 may be any double precision matrix of arbitrary size. 76 The number of columns in m2 must match the number of 77 columns in m1. 78 79 nrow1 is the number of rows in both m1 and mout. 80 81 nc1c2 i the number of columns in m1 and (by necessity) the 82 number of columns of m2. 83 84 nrow2 is the number of rows in both m2 and the number of columns 85 in mout. 86 87 -Detailed_Output 88 89 mout is the product matrix given by 90 91 t 92 mout = (m1) x (m2) 93 94 95 where the superscript "t" denotes the transpose matrix. 96 This is a double precision matrix of dimension nrow1 x 97 nrow2. 98 99 mout may overwrite m1 or m2. Note that this capability 100 does not exist in the Fortran version of SPICELIB; in the 101 Fortran version, the output must not overwrite either 102 input. 103 104 -Parameters 105 106 None. 107 108 -Particulars 109 110 The code reflects precisely the following mathematical expression 111 112 For each value of the subscript i from 1 to nrow1, and j from 1 113 to nrow2: 114 115 mout(i,j) = summation from k=1 to nc1c2 of ( m1(i,k) * m2(j,k) ) 116 117 Notice that the order of the subscripts of m2 are reversed from 118 what they would be if this routine merely multiplied m1 and m2. 119 It is this transposition of subscripts that makes this routine 120 multiply m1 and the TRANPOSE of m2. 121 122 -Examples 123 124 1) Let m1 = 125 126 | 1.0 2.0 3.0 | 127 | | 128 | 3.0 2.0 1.0 | 129 130 Let m2 = 131 132 | 1.0 2.0 0.0 | 133 | | 134 | 2.0 1.0 2.0 | 135 | | 136 | 1.0 2.0 0.0 | 137 | | 138 | 2.0 1.0 2.0 | 139 140 Here 141 142 nrow1 = 2 143 nc1c2 = 3 144 nrow2 = 4 145 146 147 so the call 148 149 mxmtg_c ( m1, m2, nrow1, nc1c2, nrow2, mout ); 150 151 152 produces the matrix 153 154 155 mout = | 5.0 10.0 5.0 10.0 | 156 | | 157 | 7.0 10.0 7.0 10.0 | 158 159 160 -Restrictions 161 162 No error checking is performed to prevent numeric overflow or 163 underflow. 164 165 No error checking is performed to determine if the input and 166 output matrices have, in fact, been correctly dimensioned. 167 168 The user is responsible for checking the magnitudes of the 169 elements of m1 and m2 so that a floating point overflow does 170 not occur. 171 172 -Exceptions 173 174 Error free. 175 176 -Files 177 178 None. 179 180 -Author_and_Institution 181 182 N.J. Bachman (JPL) 183 W.M. Owen (JPL) 184 185 -Literature_References 186 187 None. 188 189 -Version 190 191 -CSPICE Version 1.2.0, 28-AUG-2001 (NJB) 192 193 Const-qualified input arrays. 194 195 -CSPICE Version 1.1.0, 08-FEB-1998 (NJB) 196 197 Corrected a comment describing the local macro INDEX. Made 198 miscellaneous code format corrections. 199 200 -CSPICE Version 1.0.0, 25-OCT-1997 (NJB) 201 202 Based on SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) 203 204 -Index_Entries 205 206 matrix times matrix_transpose n-dimensional_case 207 208 -& 209 */ 210 211 { /* Begin mxmtg_c */ 212 213 214 215 /* 216 Local macros 217 218 We'd like to be able to refer to the elements of the input and output 219 matrices using normal subscripts, for example, m1[2][3]. Since the 220 compiler doesn't know how to compute index offsets for the array 221 arguments, which have user-adjustable size, we must compute the 222 offsets ourselves. To make syntax a little easier to read (we hope), 223 we'll use macros to do the computations. 224 225 The macro INDEX(width, i,j) computes the index offset from the array 226 base of the element at position [i][j] in a 2-dimensional matrix 227 having the number of columns indicated by width. For example, if 228 the input matrix m1 has 2 rows and 3 columns, the element at position 229 [0][1] would be indicated by 230 231 m1[ INDEX(3,0,1) ] 232 233 */ 234 235 #define INDEX( width, row, col ) ( (row)*(width) + (col) ) 236 237 238 /* 239 Local variables 240 */ 241 SpiceDouble innerProduct; 242 SpiceDouble *tmpmat; 243 SpiceDouble *loc_m1; 244 SpiceDouble *loc_m2; 245 246 SpiceInt col; 247 SpiceInt nelts; 248 249 SpiceInt row; 250 SpiceInt i; 251 252 size_t size; 253 254 255 /* 256 Allocate space for a temporary copy of the output matrix, which 257 has nrow1 rows and nc1c2 columns. 258 */ 259 nelts = nrow1 * nrow2; 260 261 size = (size_t) ( nelts * sizeof(SpiceDouble) ); 262 263 tmpmat = (SpiceDouble *) malloc ( size ); 264 265 if ( tmpmat == (SpiceDouble *)0 ) 266 { 267 chkin_c ( "mxmtg_c" ); 268 setmsg_c ( "An attempt to create a temporary matrix failed." ); 269 sigerr_c ( "SPICE(MEMALLOCFAILED)" ); 270 chkout_c ( "mxmtg_c" ); 271 return; 272 } 273 274 /* 275 Cast the input pointers to pointers to SpiceDoubles. Note: the 276 original variables are pointers to void so that callers may 277 supply the array names as arguments without casting them to 278 SpiceDoubles. The naked array name is considered by the compiler 279 to be an incompatible pointer type with (SpiceDouble *), so we 280 can't simply declare the arguments to be (SpiceDouble *). On the 281 other hand, every pointer type can be cast to (void *). 282 */ 283 284 loc_m1 = (SpiceDouble *) m1; 285 loc_m2 = (SpiceDouble *) m2; 286 287 288 /* 289 Compute the product. The matrix element at position (row,col) is 290 the inner product of the row of m1 having index row and the 291 row of m2 having index col. We compute index offsets using 292 the macro INDEX. 293 */ 294 295 for ( row = 0; row < nrow1; row++ ) 296 { 297 298 for ( col = 0; col < nrow2; col++ ) 299 { 300 innerProduct = 0.0; 301 302 for ( i = 0; i < nc1c2; i++ ) 303 { 304 innerProduct += loc_m1[ INDEX(nc1c2, row, i) ] 305 * loc_m2[ INDEX(nc1c2, col, i) ]; 306 } 307 308 tmpmat [ INDEX( nrow2, row, col ) ] = innerProduct; 309 } 310 } 311 312 /* 313 Move the result from tmpmat into mout. 314 */ 315 MOVED ( tmpmat, nelts, mout ); 316 317 /* 318 Free the temporary matrix. 319 */ 320 free ( tmpmat ); 321 322 323 } /* End mxmtg_c */ 324