1 /* 2 3 -Procedure mtxvg_c ( Matrix transpose times vector, general dimension ) 4 5 -Abstract 6 7 Multiply the transpose of a matrix and a vector of arbitrary size. 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, VECTOR 41 42 */ 43 #include <stdlib.h> 44 #include "SpiceUsr.h" 45 #include "SpiceZmc.h" 46 #include "SpiceZfc.h" 47 #include "SpiceZim.h" 48 #undef mtxvg_c 49 50 mtxvg_c(const void * m1,const void * v2,SpiceInt ncol1,SpiceInt nr1r2,void * vout)51 void mtxvg_c ( const void * m1, 52 const void * v2, 53 SpiceInt ncol1, 54 SpiceInt nr1r2, 55 void * vout ) 56 /* 57 58 -Brief_I/O 59 60 VARIABLE I/O DESCRIPTION 61 -------- --- -------------------------------------------------- 62 m1 I Left-hand matrix to be multiplied. 63 v2 I Right-hand vector to be multiplied. 64 ncol1 I Column dimension of m1 and length of vout. 65 nr1r2 I Row dimension of m1 and length of v2. 66 vout O Product vector m1 transpose * v2. 67 68 -Detailed_Input 69 70 m1 is a double precision matrix of arbitrary size which 71 forms the left-hand matrix of the multiplication. 72 73 v2 is a double precision vector on the right of the 74 multiplication. 75 76 ncol1 is the column dimension of m1 and length of vout. 77 78 nr1r2 is the row dimension of m1 and length of v2. 79 80 -Detailed_Output 81 82 vout is the double precision vector which results from 83 the multiplication 84 85 t 86 vout = (m1) x v2 87 88 where the superscript t denotes the transpose of a matrix. 89 vout has length ncol1. 90 91 vout may overwrite m1 or v2. Note that this capability 92 does not exist in the Fortran version of SPICELIB; in the 93 Fortran version, the output must not overwrite either 94 input. 95 96 -Parameters 97 98 None. 99 100 -Particulars 101 102 The code reflects precisely the following mathematical expression 103 104 For each value of the subscript i from 1 to ncol1, 105 106 vout(i) = Summation from k=1 to nr1r2 of ( m1(k,i) * v2(k) ) 107 108 -Examples 109 110 1) Suppose that 111 112 | 1 2 | 113 m1 = | 1 3 | 114 | 1 4 | 115 116 117 and that 118 119 | 1 | 120 v2 = | 2 | 121 | 3 | 122 123 124 Then calling mxvg_c as shown 125 126 mtxvg_c ( m1, v2, 2, 3, vout ); 127 128 129 will yield the following vector value for vout: 130 131 vout = | 6 | 132 | 20 | 133 134 -Restrictions 135 136 1) The user is responsible for checking the magnitudes of the 137 elements of m1 and v2 so that a floating point overflow does 138 not occur. 139 140 -Exceptions 141 142 Error free. 143 144 -Files 145 146 None. 147 148 -Author_and_Institution 149 150 N.J. Bachman (JPL) 151 W.M. Owen (JPL) 152 153 -Literature_References 154 155 None. 156 157 -Version 158 159 -CSPICE Version 1.2.0, 28-AUG-2001 (NJB) 160 161 Const-qualified input arrays. 162 163 -CSPICE Version 1.1.0, 08-FEB-1998 (NJB) 164 165 Corrected a comment describing the local macro INDEX. Made 166 miscellaneous code format corrections. 167 168 Based on SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) 169 170 -Index_Entries 171 172 matrix transpose times n-dimensional vector 173 174 -& 175 */ 176 177 { /* Begin mxvg_c */ 178 179 /* 180 Local macros 181 182 We'd like to be able to refer to the elements of the input and output 183 matrices using normal subscripts, for example, m1[2][3]. Since the 184 compiler doesn't know how to compute index offsets for the array 185 arguments, which have user-adjustable size, we must compute the 186 offsets ourselves. To make syntax a little easier to read (we hope), 187 we'll use macros to do the computations. 188 189 The macro INDEX(width, i,j) computes the index offset from the array 190 base of the element at position [i][j] in a 2-dimensional matrix 191 having the number of columns indicated by width. For example, if 192 the input matrix m1 has 2 rows and 3 columns, the element at position 193 [0][1] would be indicated by 194 195 m1[ INDEX(3,0,1) ] 196 197 */ 198 199 #define INDEX( width, row, col ) ( (row)*(width) + (col) ) 200 201 202 /* 203 Local variables 204 */ 205 SpiceDouble innerProduct; 206 SpiceDouble *tmpvec; 207 SpiceDouble *loc_m1; 208 SpiceDouble *loc_v2; 209 210 SpiceInt row; 211 SpiceInt i; 212 213 size_t size; 214 215 216 /* 217 Allocate space for a temporary copy of the output vector, which 218 has ncol1 rows. 219 */ 220 size = (size_t) ( ncol1 * sizeof(SpiceDouble) ); 221 222 tmpvec = (SpiceDouble *) malloc ( size ); 223 224 if ( tmpvec == (SpiceDouble *)0 ) 225 { 226 chkin_c ( "mtxvg_c" ); 227 setmsg_c ( "An attempt to create a temporary vector failed." ); 228 sigerr_c ( "SPICE(MEMALLOCFAILED)" ); 229 chkout_c ( "mtxvg_c" ); 230 return; 231 } 232 233 /* 234 Cast the input pointers to pointers to SpiceDoubles. Note: the 235 original variables are pointers to void so that callers may 236 supply the array names as arguments without casting them to 237 SpiceDoubles. The naked array name is considered by the compiler 238 to be an incompatible pointer type with (SpiceDouble *), so we 239 can't simply declare the arguments to be (SpiceDouble *). On the 240 other hand, every pointer type can be cast to (void *). 241 */ 242 243 loc_m1 = (SpiceDouble *) m1; 244 loc_v2 = (SpiceDouble *) v2; 245 246 247 /* 248 Compute the product. The vector element at position (row) is 249 the inner product of the column of m1 having index row and v2. 250 We compute index offsets using the macro INDEX. 251 */ 252 253 for ( row = 0; row < ncol1; row++ ) 254 { 255 256 innerProduct = 0.0; 257 258 for ( i = 0; i < nr1r2; i++ ) 259 { 260 innerProduct += loc_m1[ INDEX(ncol1, i, row ) ] * loc_v2[i]; 261 } 262 263 tmpvec [ row ] = innerProduct; 264 } 265 266 /* 267 Move the result from tmpvec into vout. 268 */ 269 MOVED ( tmpvec, ncol1, vout ); 270 271 /* 272 Free the temporary vector. 273 */ 274 free ( tmpvec ); 275 276 277 } /* End mtxvg_c */ 278