00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include <stdio.h>
00045 #include "sp_matrix.h"
00046
00052
00053 MATRIX create_matrix (int rows, int cols)
00054
00055
00056
00057
00058 {
00059 MATRIX m;
00060
00061 MROWS (m) = rows;
00062 MCOLS (m) = cols;
00063
00064 {
00065 int i, j;
00066
00067 for (i = 0; i < MROWS (m); i++)
00068 for (j = 0; j < MCOLS (m); j++)
00069 MDATA (m, i, j) = 0;
00070 }
00071
00072 return m;
00073 }
00074
00075
00076 void initialize_matrix (MATRIX *m, int rows, int cols)
00077
00078
00079
00080
00081 {
00082 MROWS (*m) = rows;
00083 MCOLS (*m) = cols;
00084
00085 {
00086 int i, j;
00087
00088 for (i = 0; i < MROWS (*m); i++)
00089 for (j = 0; j < MCOLS (*m); j++)
00090 MDATA (*m, i, j) = 0;
00091 }
00092
00093 }
00094
00095
00096
00097 void print_matrix (char *message, MATRIX const *m)
00098
00099
00100
00101
00102 {
00103 int i, j;
00104
00105 printf ("%s\n",message);
00106 printf("%d %d \n",MROWS (*m),MCOLS (*m));
00107 if ((MROWS (*m) <= MAX_ROWS) && (MCOLS (*m) <= MAX_COLS))
00108 for (i = 0; i < MROWS (*m); i++)
00109 {
00110 for (j = 0; j < MCOLS (*m); j++)
00111 printf ("%10.5f ", MDATA (*m, i, j));
00112 printf ("\n");
00113 }
00114 else printf ("Dimension incorrecta!");
00115 printf ("\n");
00116 }
00117
00118
00119 VECTOR create_vector (int elements)
00120
00121
00122
00123
00124 {
00125 VECTOR v;
00126
00127 VELEMENTS (v) = elements;
00128
00129 {
00130 int i;
00131
00132 for (i = 0; i < VELEMENTS (v); i++)
00133 VDATA (v, i) = 0;
00134 }
00135
00136 return v;
00137 }
00138
00139
00140 void initialize_vector (VECTOR *v, int elements)
00141
00142
00143
00144
00145 {
00146 VELEMENTS (*v) = elements;
00147
00148 {
00149 int i;
00150
00151 for (i = 0; i < VELEMENTS (*v); i++)
00152 VDATA (*v, i) = 0;
00153 }
00154 }
00155
00156
00157 void print_vector (char *message, VECTOR const *v)
00158
00159
00160
00161
00162 {
00163 int i;
00164
00165 printf ("%s\n",message);
00166 if (VELEMENTS (*v) <= MAX_ROWS)
00167 for (i = 0; i < VELEMENTS (*v); i++)
00168 {
00169 printf ("%f ", VDATA (*v, i));
00170 printf ("\n");
00171 }
00172 else printf ("Dimension incorrecta!");
00173 printf ("\n");
00174 }
00175
00176
00177 float cross_product (MATRIX const *m, int f1, int c1, int f2, int c2)
00178
00179
00180
00181 {
00182 return MDATA (*m, f1, c1) * MDATA (*m, f2, c2) - MDATA (*m, f1, c2) * MDATA (*m, f2, c1);
00183 }
00184
00185
00186 int determinant (MATRIX const *m, float *result)
00187
00188
00189 {
00190 if (!M_SQUARE (*m))
00191 {
00192 printf ("ERROR (determinant): MATRIX must be square!\n");
00193 print_matrix ("MATRIX:", m);
00194 return -1;
00195 }
00196 else
00197 {
00198
00199 if (MROWS (*m) == 1)
00200 *result = MDATA (*m, 0, 0);
00201 else if (MROWS (*m) == 2)
00202 *result = cross_product (m, 0, 0, 1, 1);
00203 else
00204 *result = MDATA (*m, 0, 0) * cross_product (m, 1, 1, 2, 2)
00205 - MDATA (*m, 0, 1) * cross_product (m, 1, 0, 2, 2)
00206 + MDATA (*m, 0, 2) * cross_product (m, 1, 0, 2, 1);
00207
00208
00209 return 1;
00210 }
00211 }
00212
00213
00214 int inverse_matrix (MATRIX const *m, MATRIX *n)
00215
00216
00217
00218 {
00219 if (!M_SQUARE (*m))
00220 {
00221 printf ("ERROR (inverse_matrix): MATRIX must be square!\n");
00222 print_matrix ("MATRIX:", m);
00223 n->cols=0; n->rows=0;
00224 return -1;
00225 }
00226 else
00227 {
00228 float det;
00229 int res;
00230
00231 res = determinant (m,&det);
00232
00233 if (res == -1)
00234 {
00235 printf ("ERROR (inverse_matrix): singular MATRIX!\n");
00236 print_matrix ("MATRIX:", m);
00237 return -1;
00238 }
00239 else
00240 {
00241 initialize_matrix (n, MROWS (*m), MCOLS (*m));
00242 if (MROWS (*m) == 1)
00243 {
00244 MDATA (*n, 0, 0) = 1 / det ;
00245 }
00246 else if (MROWS (*m) == 2)
00247 {
00248 MDATA (*n, 0, 0) = MDATA (*m, 1, 1) / det ;
00249 MDATA (*n, 0, 1) = -MDATA (*m, 0, 1) / det ;
00250 MDATA (*n, 1, 0) = -MDATA (*m, 1, 0) / det ;
00251 MDATA (*n, 1, 1) = MDATA (*m, 0, 0) / det ;
00252 }
00253 else
00254 {
00255 MDATA (*n, 0, 0) = cross_product (m, 1, 1, 2, 2) / det ;
00256 MDATA (*n, 0, 1) = -cross_product (m, 0, 1, 2, 2) / det ;
00257 MDATA (*n, 0, 2) = cross_product (m, 0, 1, 1, 2) / det ;
00258 MDATA (*n, 1, 0) = -cross_product (m, 1, 0, 2, 2) / det ;
00259 MDATA (*n, 1, 1) = cross_product (m, 0, 0, 2, 2) / det ;
00260 MDATA (*n, 1, 2) = -cross_product (m, 0, 0, 1, 2) / det ;
00261 MDATA (*n, 2, 0) = cross_product (m, 1, 0, 2, 1) / det ;
00262 MDATA (*n, 2, 1) = -cross_product (m, 0, 0, 2, 1) / det ;
00263 MDATA (*n, 2, 2) = cross_product (m, 0, 0, 1, 1) / det ;
00264 }
00265 }
00266 }
00267 return 1;
00268 }
00269
00270
00271 int multiply_matrix_vector (MATRIX const *m, VECTOR const *v, VECTOR *r)
00272
00273
00274
00275
00276 {
00277 if (! (MV_COMPAT_DIM (*m, *v)))
00278 {
00279 printf ("ERROR (multiply_matrix_vector): MATRIX and VECTOR dimensions incompatible!\n");
00280 print_matrix ("MATRIX:", m);
00281 print_vector ("VECTOR:", v);
00282 return -1;
00283 }
00284 else
00285 {
00286 int i, j;
00287 float datum;
00288
00289 VELEMENTS (*r) = MROWS (*m);
00290
00291 for (i = 0; i < MROWS (*m); i++)
00292 {
00293 datum = 0;
00294 for (j = 0; j < VELEMENTS (*v); j++)
00295 datum = datum + MDATA (*m, i, j) * VDATA (*v, j);
00296 VDATA (*r, i) = datum;
00297 }
00298 }
00299 return 1;
00300 }
00301