1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libann/matrix.c,v 1.9 2011/02/23 13:02:19 mattaboni Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2008
7 *
8 * Pierangelo Masarati <masarati@aero.polimi.it>
9 * Paolo Mantegazza <mantegazza@aero.polimi.it>
10 *
11 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12 * via La Masa, 34 - 20156 Milano, Italy
13 * http://www.aero.polimi.it
14 *
15 * Changing this copyright notice is forbidden.
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation (version 2 of the License).
20 *
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30 */
31 /*
32 * Copyright (C) 2010
33 *
34 * Mattia Mattaboni <mattaboni@aero.polimi.it>
35 */
36
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
40 #include "matrix.h"
41
42 /* inizializza un elemento della classe matrice */
matrix_init(matrix * MAT,unsigned Nrow,unsigned Ncolumn)43 mat_res_t matrix_init( matrix *MAT, unsigned Nrow, unsigned Ncolumn ){
44
45 unsigned i;
46
47 if( Nrow<= 0 || Ncolumn <= 0 ){
48 matrix_error( MAT_DIMENSION, "matrix_init" );
49 return MAT_DIMENSION;
50 }
51
52 MAT->Nrow = Nrow;
53 MAT->Ncolumn = Ncolumn;
54 if( !( MAT->mat = (double **)calloc( Nrow, sizeof(double *) ) ) ){
55 matrix_error( MAT_NO_MEMORY, "matrix_init" );
56 return MAT_NO_MEMORY;
57 }
58 for( i=0; i<Nrow; i++ ){
59 if( !( MAT->mat[i] = (double *)calloc( Ncolumn, sizeof(double) ) ) ){
60 matrix_error( MAT_NO_MEMORY, "matrix_init" );
61 return MAT_NO_MEMORY;
62 }
63 }
64
65 return MAT_OK;
66 }
67 /* inizializza un elemento della classe vettore */
vector_init(vector * VEC,unsigned dimension)68 mat_res_t vector_init( vector *VEC, unsigned dimension ){
69 if( dimension<= 0 ){
70 matrix_error( MAT_DIMENSION, "vector_init" );
71 return MAT_DIMENSION;
72 }
73
74 VEC->dimension = dimension;
75
76 if( !( VEC->vec = (double *)calloc( dimension, sizeof(double) ) ) ){
77 matrix_error( MAT_NO_MEMORY, "vector_init");
78 return MAT_NO_MEMORY;
79 }
80
81 return MAT_OK;
82 }
83 /* distrugge un elemento della classe matrice*/
matrix_destroy(matrix * MAT)84 mat_res_t matrix_destroy( matrix *MAT ) {
85
86 unsigned i;
87
88 for( i=0; i<MAT->Nrow; i++ ){
89 free(MAT->mat[i]);
90 }
91 free(MAT->mat);
92
93 return MAT_OK;
94 }
95
96 /* distrugge un elemento della classe vettore */
vector_destroy(vector * VEC)97 mat_res_t vector_destroy( vector *VEC ) {
98
99 free(VEC->vec);
100
101 return MAT_OK;
102 }
103
104 /* OPERAZIONI TRA MATRICI E VETTORI */
105
106 /* azzera una matrice MAT = zeros*/
matrix_null(matrix * MAT)107 mat_res_t matrix_null( matrix *MAT ){
108
109 unsigned i;
110
111 for( i=0; i<MAT->Nrow; i++ ){
112 if( !(memset( MAT->mat[i], 0, MAT->Ncolumn*sizeof(double) ) ) ){
113 matrix_error( MAT_GEN_ERROR, "matrix_null" );
114 return MAT_GEN_ERROR;
115 }
116 }
117
118 return MAT_OK;
119 }
120
121 /* Identity matrix: MAT = K*eye */
matrix_eye(matrix * MAT,double K)122 mat_res_t matrix_eye( matrix *MAT, double K ){
123
124 unsigned i;
125
126 /* azzero la matrice del risultato */
127 if( matrix_null(MAT) != MAT_OK ){
128 matrix_error( MAT_GEN_ERROR, "matrix_eye" );
129 return MAT_GEN_ERROR;
130 }
131
132 for( i=0; i<MAT->Nrow; i++ ){
133 MAT->mat[i][i] = K;
134 }
135
136 return MAT_OK;
137 }
138
139 /* copia una matrice MAT1 = MAT2*/
matrix_copy(matrix * MAT1,matrix * MAT2,double K)140 mat_res_t matrix_copy( matrix *MAT1, matrix *MAT2, double K ){
141
142 unsigned i, j;
143
144 /* controllo dimensionale */
145 if( MAT1->Ncolumn != MAT2->Ncolumn || MAT1->Nrow != MAT2->Nrow){
146 matrix_error( MAT_DIMENSION, "matrix_copy" );
147 return MAT_DIMENSION;
148 }
149
150 for( i=0; i<MAT1->Nrow; i++ ){
151 for( j=0; j<MAT1->Ncolumn; j++ ){
152 MAT1->mat[i][j] = K*MAT2->mat[i][j];
153 }
154 }
155
156 return MAT_OK;
157 }
158 /* copia un vettore VEC1 = VEC2*/
vector_copy(vector * VEC1,vector * VEC2,double K)159 mat_res_t vector_copy( vector *VEC1, vector *VEC2, double K ){
160
161 unsigned i;
162
163 /* controllo dimensionale */
164 if( VEC1->dimension != VEC2->dimension){
165 matrix_error( MAT_DIMENSION, "vector_copy" );
166 return MAT_DIMENSION;
167 }
168
169 for( i=0; i<VEC1->dimension; i++ ){
170 VEC1->vec[i] = K*VEC2->vec[i];
171 }
172
173 return MAT_OK;
174 }
175 /* azzera un vettore VEC = zeros */
vector_null(vector * VEC)176 mat_res_t vector_null( vector *VEC ){
177
178 if( !(memset( VEC->vec, 0, VEC->dimension*sizeof(double) ) ) ){
179 matrix_error( MAT_GEN_ERROR, "vector_null" );
180 return MAT_GEN_ERROR;
181 }
182
183 return MAT_OK;
184 }
185
186 /* prodotto tra matrici MAT_R = K*MAT1*MAT2*/
matrix_prod(matrix * MAT1,matrix * MAT2,matrix * MAT_R,double K)187 mat_res_t matrix_prod( matrix *MAT1 ,matrix *MAT2, matrix *MAT_R, double K ){
188
189 unsigned i,j,k;
190
191 /* controllo dimensionale */
192 if( MAT1->Ncolumn != MAT2->Nrow || MAT_R->Nrow != MAT1->Nrow || MAT_R->Ncolumn != MAT2->Ncolumn ){
193 matrix_error( MAT_DIMENSION, "matrix_prod" );
194 return MAT_DIMENSION;
195 }
196 /* azzero la matrice del risultato */
197 if( matrix_null(MAT_R) != MAT_OK ){
198 matrix_error( MAT_GEN_ERROR, "matrix_prod" );
199 return MAT_GEN_ERROR;
200 }
201 for( i=0; i<MAT1->Nrow; i++ ){
202 for( j=0; j<MAT2->Ncolumn; j++ ){
203 for( k=0; k<MAT1->Ncolumn; k++ ){
204 MAT_R->mat[i][j] += K*MAT1->mat[i][k]*MAT2->mat[k][j];
205 }
206 }
207 }
208
209 return MAT_OK;
210 }
211 /* prodotto tra matrici MAT_R = K*MAT1*MAT2*/
matrix_prod_sym(matrix * MAT1,matrix * MAT2,matrix * MAT_R,double K)212 mat_res_t matrix_prod_sym( matrix *MAT1 ,matrix *MAT2, matrix *MAT_R, double K ){
213
214 unsigned i,j,k;
215
216 /* controllo dimensionale */
217 if( MAT1->Ncolumn != MAT2->Nrow || MAT_R->Nrow != MAT1->Nrow || MAT_R->Ncolumn != MAT2->Ncolumn ){
218 matrix_error( MAT_DIMENSION, "matrix_prod" );
219 return MAT_DIMENSION;
220 }
221 /* azzero la matrice del risultato */
222 if( matrix_null(MAT_R) != MAT_OK ){
223 matrix_error( MAT_GEN_ERROR, "matrix_prod" );
224 return MAT_GEN_ERROR;
225 }
226 for( i=0; i<MAT1->Nrow; i++ ){
227 for( j=i; j<MAT2->Ncolumn; j++ ){
228 for( k=0; k<MAT1->Ncolumn; k++ ){
229 MAT_R->mat[i][j] += K*MAT1->mat[i][k]*MAT2->mat[k][j];
230 }
231 MAT_R->mat[j][i] = MAT_R->mat[i][j];
232 }
233 }
234
235 return MAT_OK;
236 }
237
238 /* matrice trasposta MAT1 = MAT2^T */
matrix_transpose(matrix * MAT1,matrix * MAT2)239 mat_res_t matrix_transpose( matrix *MAT1 ,matrix *MAT2){
240
241 unsigned i,j;
242
243 /* controllo dimensionale */
244 if( MAT1->Ncolumn != MAT2->Nrow || MAT1->Nrow != MAT2->Ncolumn ){
245 matrix_error( MAT_DIMENSION, "matrix_transpose" );
246 return MAT_DIMENSION;
247 }
248
249 for( i=0; i<MAT2->Nrow; i++ ){
250 for( j=0; j<MAT2->Ncolumn; j++ ){
251 MAT1->mat[j][i] = MAT2->mat[i][j];
252 }
253 }
254
255 return MAT_OK;
256 }
257 /* prodotto tra matrice e matrice trasposta MAT_R =K*MAT1^T*MAT2*/
matrix_transpose_prod(matrix * MAT1,matrix * MAT2,matrix * MAT_R,double K)258 mat_res_t matrix_transpose_prod( matrix *MAT1 ,matrix *MAT2, matrix *MAT_R, double K ){
259
260 unsigned i,j,k;
261
262 /* controllo dimensionale */
263 if( MAT1->Nrow != MAT2->Nrow || MAT_R->Nrow != MAT1->Ncolumn || MAT_R->Ncolumn != MAT2->Ncolumn ){
264 matrix_error( MAT_DIMENSION, "matrix_transpose_prod" );
265 return MAT_DIMENSION;
266 }
267 /* azzero la matrice del risultato */
268 if( matrix_null(MAT_R) != MAT_OK ){
269 matrix_error( MAT_GEN_ERROR, "matrix_transpose_prod" );
270 return MAT_GEN_ERROR;
271 }
272 for( i=0; i<MAT1->Ncolumn; i++ ){
273 for( j=0; j<MAT2->Ncolumn; j++ ){
274 for( k=0; k<MAT1->Nrow; k++ ){
275 MAT_R->mat[i][j] += K*MAT1->mat[k][i]*MAT2->mat[k][j];
276 }
277 }
278 }
279
280 return MAT_OK;
281 }
282
283 /* prodotto tra matrice e matrice trasposta MAT_R =K*MAT1*MAT2^T*/
matrix_prod_transpose(matrix * MAT1,matrix * MAT2,matrix * MAT_R,double K)284 mat_res_t matrix_prod_transpose( matrix *MAT1 ,matrix *MAT2, matrix *MAT_R, double K ){
285
286 unsigned i,j,k;
287
288 /* controllo dimensionale */
289 if( MAT1->Ncolumn != MAT2->Ncolumn || MAT_R->Nrow != MAT1->Nrow || MAT_R->Ncolumn != MAT2->Nrow ){
290 matrix_error( MAT_DIMENSION, "matrix_prod_transpose" );
291 return MAT_DIMENSION;
292 }
293 /* azzero la matrice del risultato */
294 if( matrix_null(MAT_R) != MAT_OK ){
295 matrix_error( MAT_GEN_ERROR, "matrix_prod_transpose" );
296 return MAT_GEN_ERROR;
297 }
298 for( i=0; i<MAT1->Nrow; i++ ){
299 for( j=0; j<MAT2->Nrow; j++ ){
300 for( k=0; k<MAT1->Ncolumn; k++ ){
301 MAT_R->mat[i][j] += K*MAT1->mat[i][k]*MAT2->mat[j][k];
302 }
303 }
304 }
305
306 return MAT_OK;
307 }
308
309 /* prodotto scalare RES = VEC1 dot VEC2*/
scalar_prod(vector * VEC1,vector * VEC2,double * RES)310 mat_res_t scalar_prod( vector *VEC1, vector *VEC2, double *RES){
311
312 unsigned i;
313 double res;
314
315 /* controllodimensionale */
316 if( VEC1->dimension != VEC2->dimension ){
317 matrix_error( MAT_DIMENSION, "matrix_vector_prod" );
318 return MAT_DIMENSION;
319 }
320
321 res = 0.;
322 for( i=0; i<VEC1->dimension; i++){
323 res += VEC1->vec[i]*VEC2->vec[i];
324 }
325
326 *RES = res;
327
328 return MAT_OK;
329
330 }
331
332 /* prodotto vettore*vettore = matrice RES = K*VEC1*VEC2^T */
vector_vector_prod(vector * VEC1,vector * VEC2,matrix * RES,double K)333 mat_res_t vector_vector_prod( vector *VEC1, vector *VEC2, matrix *RES, double K){
334
335 unsigned i, j;
336
337 /* controllodimensionale */
338 if( VEC1->dimension != RES->Nrow || VEC2->dimension != RES->Ncolumn){
339 matrix_error( MAT_DIMENSION, "vector_vector_prod" );
340 return MAT_DIMENSION;
341 }
342
343 for( i=0; i<VEC1->dimension; i++){
344 for( j=0; j<VEC2->dimension; j++){
345 RES->mat[i][j] = K*VEC1->vec[i]*VEC2->vec[j];
346 }
347 }
348
349 return MAT_OK;
350
351 }
352
353 /* prodotto matrice vettore VEC_R = MAT*VEC */
matrix_vector_prod(matrix * MAT,vector * VEC,vector * VEC_R)354 mat_res_t matrix_vector_prod( matrix *MAT, vector *VEC, vector *VEC_R){
355
356 unsigned i,j;
357
358 /* controllodimensionale */
359 if( MAT->Ncolumn != VEC->dimension || MAT->Nrow != VEC_R->dimension ){
360 matrix_error( MAT_DIMENSION, "matrix_vector_prod" );
361 return MAT_DIMENSION;
362 }
363 /* azzero il vettore risultato */
364 if( vector_null(VEC_R) != MAT_OK ){
365 matrix_error( MAT_GEN_ERROR, "matrix_vector_prod");
366 return MAT_GEN_ERROR;
367 }
368
369 for( i=0; i<MAT->Nrow; i++ ){
370 for( j=0; j<MAT->Ncolumn; j++ ){
371 VEC_R->vec[i] += MAT->mat[i][j]*VEC->vec[j];
372 }
373 }
374
375 return MAT_OK;
376 }
377
378 /* prodotto tra una matrice trasposta ed un vettore VEC_R = MAT^T*VEC */
matrixT_vector_prod(matrix * MAT,vector * VEC,vector * VEC_R)379 mat_res_t matrixT_vector_prod( matrix *MAT, vector *VEC, vector *VEC_R){
380
381 unsigned i,j;
382
383 /* controllo dimensionale */
384 if( MAT->Nrow != VEC->dimension || MAT->Ncolumn != VEC_R->dimension ){
385 matrix_error( MAT_DIMENSION, "matrixT_vector_prod" );
386 return MAT_DIMENSION;
387 }
388 /* annullo il vettore risulatante */
389 if( vector_null(VEC_R) != MAT_OK ){
390 matrix_error( MAT_GEN_ERROR, "matrixT_vector_prod" );
391 return MAT_GEN_ERROR;
392 }
393
394 //for( i=0; i<MAT->Nrow; i++ ){
395 // for( j=0; j<MAT->Ncolumn; j++ ){
396 // VEC_R->vec[j] += MAT->mat[i][j]*VEC->vec[i];
397 // }
398 //}
399 for( i=0; i<MAT->Ncolumn; i++ ){
400 for( j=0; j<MAT->Nrow; j++ ){
401 VEC_R->vec[i] += MAT->mat[j][i]*VEC->vec[j];
402 }
403 }
404
405 return MAT_OK;
406 }
407
408 /* somma tra matrici MAT_R = MAT1+K*MAT */
matrix_sum(matrix * MAT1,matrix * MAT2,matrix * MAT_R,double K)409 mat_res_t matrix_sum( matrix *MAT1, matrix *MAT2, matrix *MAT_R, double K ){
410
411 unsigned i,j;
412
413 /* controllo dimensionale */
414 if( MAT1->Ncolumn != MAT2->Ncolumn || MAT_R->Ncolumn != MAT1->Ncolumn || MAT_R->Nrow != MAT2->Nrow || MAT1->Nrow != MAT2->Nrow ){
415 matrix_error( MAT_DIMENSION, "matrix_sum" );
416 return MAT_DIMENSION;
417 }
418 for( i=0; i<MAT1->Nrow; i++ ){
419 for( j=0; j<MAT1->Ncolumn; j++ ){
420 MAT_R->mat[i][j] = MAT1->mat[i][j] + K*MAT2->mat[i][j];
421 }
422 }
423
424 return MAT_OK;
425 }
426
427 /* somma tra matrici MAT_R = MAT1+K*MAT2^T */
matrix_sum_transpose(matrix * MAT1,matrix * MAT2,matrix * MAT_R,double K)428 mat_res_t matrix_sum_transpose( matrix *MAT1, matrix *MAT2, matrix *MAT_R, double K ){
429
430 unsigned i,j;
431
432 /* controllo dimensionale */
433 if( MAT1->Ncolumn != MAT2->Nrow || MAT_R->Ncolumn != MAT1->Ncolumn || MAT_R->Nrow != MAT1->Nrow || MAT1->Nrow != MAT2->Ncolumn ){
434 matrix_error( MAT_DIMENSION, "matrix_sum_transpose" );
435 return MAT_DIMENSION;
436 }
437 for( i=0; i<MAT1->Nrow; i++ ){
438 for( j=0; j<MAT1->Ncolumn; j++ ){
439 MAT_R->mat[i][j] = MAT1->mat[i][j] + K*MAT2->mat[j][i];
440 }
441 }
442
443 return MAT_OK;
444 }
445
446 /* somma tra vettori VEC_R = VEC1+K*VEC2 */
vector_sum(vector * VEC1,vector * VEC2,vector * VEC_R,double K)447 mat_res_t vector_sum( vector *VEC1, vector *VEC2, vector *VEC_R, double K ){
448
449 unsigned i;
450
451 /* controllo dimensionale */
452 if( VEC1->dimension != VEC2->dimension || VEC_R->dimension != VEC1->dimension ){
453 matrix_error( MAT_DIMENSION, "vector_sum" );
454 return MAT_DIMENSION;
455 }
456 for( i=0; i<VEC1->dimension; i++ ){
457 VEC_R->vec[i] = VEC1->vec[i] + K*VEC2->vec[i];
458 }
459
460 return MAT_OK;
461 }
462
463
464 /* FUNZIONI ACCESSORIE */
465
466 /* scrive a video o su file una matrice */
matrix_write(matrix * MAT,FILE * fh,unsigned flags)467 mat_res_t matrix_write( matrix *MAT, FILE *fh, unsigned flags ){
468
469 unsigned i,j;
470
471 if( flags & W_M_TEXT ) fprintf( fh, "matrix = [\n" );
472 if( flags & W_M_BIN ) fprintf( fh, "\n" );
473 for( i=0; i<MAT->Nrow; i++ ){
474 for( j=0; j<MAT->Ncolumn; j++ ){
475 fprintf( fh, "%15.16e ", MAT->mat[i][j] );
476 }
477 fprintf( fh, "\n");
478 }
479
480 if( flags & W_M_TEXT ) fprintf( fh, "]\n" );
481 if( flags & W_M_BIN ) fprintf( fh, "\n" );
482
483 return MAT_OK;
484 }
485
486 /* scrive a video o su file un vettore */
vector_write(vector * VEC,FILE * fh,unsigned flags)487 mat_res_t vector_write( vector *VEC, FILE *fh, unsigned flags ){
488
489 unsigned i;
490
491 if( flags & W_M_TEXT ) fprintf( fh, "vector = [\n" );
492 /*if( flags & W_M_BIN ) fprintf( fh, "\n" );*/
493 for( i=0; i<VEC->dimension; i++ ){
494 if( flags & W_M_BIN ) fprintf( fh, "\n" );
495 fprintf( fh, "%15.16e ", VEC->vec[i] );
496 }
497
498 if( flags & W_M_TEXT ) fprintf( fh, "]\n" );
499 if( flags & W_M_BIN ) fprintf( fh, "\n" );
500 if( flags & W_M_BIN_ROW ) fprintf( fh, "\n" );
501
502 return MAT_OK;
503 }
504
505 /* legge una matrice da file */
matrix_read(matrix * MAT,FILE * fh,unsigned flags)506 mat_res_t matrix_read( matrix *MAT, FILE * fh, unsigned flags ){
507
508 unsigned i,j;
509
510 for( i=0; i<MAT->Nrow; i++ ){
511 for( j=0; j<MAT->Ncolumn; j++ ){
512 if( fscanf( fh, "%le", &MAT->mat[i][j] ) <= 0 ){
513 return MAT_INPUT;
514 }
515 }
516 }
517
518 return MAT_OK;
519 }
520
521 /* legge una vettore da file */
vector_read(vector * VEC,FILE * fh,unsigned flags)522 mat_res_t vector_read( vector *VEC, FILE * fh, unsigned flags ){
523
524 unsigned i;
525
526 for( i=0; i<VEC->dimension; i++ ){
527 if( fscanf( fh, "%le", &VEC->vec[i] ) <= 0 ){
528 return MAT_INPUT;
529 }
530 }
531
532 return MAT_OK;
533 }
534 /* gestione degli errori */
matrix_error(mat_res_t error,const char * string)535 void matrix_error( mat_res_t error, const char * string){
536
537 switch(error) {
538 case MAT_NO_MEMORY: fprintf( stderr, "Memory error( @ %s )\n", string );
539 break;
540 case MAT_DIMENSION: fprintf( stderr, "Matrix dimension mismatch( @ %s )\n", string );
541 break;
542 case MAT_INPUT: fprintf( stderr, "Reading error( @ %s )\n",string );
543 break;
544 case MAT_GEN_ERROR: fprintf( stderr, "Error( @ %s )\n",string );
545 break;
546 default: break;
547 }
548 }
549
550 /* genera una matrice random di numeri compresi tra min e max*/
551
matrix_random(matrix * MAT,double min,double max)552 mat_res_t matrix_random( matrix *MAT, double min, double max ){
553
554 double y;
555 unsigned i,j;
556
557 for( i=0; i<MAT->Nrow; i++ ){
558 for( j=0; j<MAT->Ncolumn; j++ ){
559 y = rand();
560 y = y/RAND_MAX;
561 y = y*(max-min);
562 y += min ;
563 MAT->mat[i][j] = y;
564 }
565 }
566
567 return MAT_OK;
568 }
569
vector_random(vector * VEC,double min,double max)570 mat_res_t vector_random( vector *VEC, double min, double max ){
571
572 double y;
573 unsigned i;
574
575 for( i=0; i<VEC->dimension; i++ ){
576 y = rand();
577 y = y/RAND_MAX;
578 y = y*(max-min);
579 y += min ;
580 VEC->vec[i] = y;
581 }
582 return MAT_OK;
583 }
584
585 /* calcola il valor medio della colonna "column" della matrice "MAT" */
mean_value(matrix * MAT,int column)586 double mean_value( matrix *MAT, int column ){
587
588 unsigned i;
589 double mean = 0.;
590
591 for( i=0; i<MAT->Nrow; i++ ){
592 mean += MAT->mat[i][column];
593 }
594 return( mean/MAT->Nrow );
595 }
596
597 /* calcola la varianza della colonna "column" della matrice "MAT" */
variance(matrix * MAT,int column)598 double variance( matrix *MAT, int column ){
599
600 unsigned i;
601 double mean, var = 0.;
602
603 mean = mean_value( MAT, column);
604 for( i=0; i<MAT->Nrow; i++ ){
605 var += (MAT->mat[i][column]-mean)*(MAT->mat[i][column]-mean);
606 }
607
608 return( var/MAT->Nrow );
609 }
610
611 /* calcola il valor massimo della colonna "column" della matrice "MAT" */
maximum(matrix * MAT,int column)612 double maximum( matrix *MAT, int column ){
613
614 unsigned i;
615 double MAX;
616
617 MAX = MAT->mat[0][column];
618 for( i=0; i<MAT->Nrow; i++ ){
619 if( MAT->mat[i][column] > MAX ){
620 MAX = MAT->mat[i][column];
621 }
622 }
623
624 return MAX;
625 }
626
627 /* calcola il valor minimo della colonna "column" della matrice "MAT" */
minimum(matrix * MAT,int column)628 double minimum( matrix *MAT, int column ){
629
630 unsigned i;
631 double MIN;
632
633 MIN = MAT->mat[0][column];
634 for( i=0; i<MAT->Nrow; i++ ){
635 if( MAT->mat[i][column] < MIN ){
636 MIN = MAT->mat[i][column];
637 }
638 }
639
640 return MIN;
641 }
642
643 /* calcola la traccia di una matrice */
matrix_trace(matrix * MAT)644 double matrix_trace( matrix *MAT ){
645
646 unsigned i;
647 double trace;
648
649 trace = 0.;
650 for( i=0; i<MAT->Nrow; i++ ){
651 trace += MAT->mat[i][i];
652 }
653
654 return trace;
655 }
sub_matrix_extract(matrix * BIG,matrix * SUB,unsigned RowIndex,unsigned ColumnIndex)656 mat_res_t sub_matrix_extract( matrix *BIG, matrix *SUB, unsigned RowIndex, unsigned ColumnIndex){
657
658 unsigned i, j;
659
660 if( (RowIndex+SUB->Nrow > BIG->Nrow) || (ColumnIndex+SUB->Ncolumn > BIG->Ncolumn) ){
661 matrix_error( MAT_DIMENSION, "sub_matrix_extract" );
662 return MAT_DIMENSION;
663 }
664
665 for( i=0; i<SUB->Nrow; i++){
666 for( j=0; j<SUB->Ncolumn; j++){
667 SUB->mat[i][j] = BIG->mat[RowIndex+i][ColumnIndex+j];
668 }
669 }
670
671 return MAT_OK;
672 }
673
sub_matrix_insert(matrix * BIG,matrix * SUB,unsigned RowIndex,unsigned ColumnIndex)674 mat_res_t sub_matrix_insert( matrix *BIG, matrix *SUB, unsigned RowIndex, unsigned ColumnIndex){
675
676 unsigned i, j;
677
678 if( (RowIndex+SUB->Nrow > BIG->Nrow) || (ColumnIndex+SUB->Ncolumn > BIG->Ncolumn) ){
679 matrix_error( MAT_DIMENSION, "sub_matrix_insert" );
680 return MAT_DIMENSION;
681 }
682
683 for( i=0; i<SUB->Nrow; i++){
684 for( j=0; j<SUB->Ncolumn; j++){
685 BIG->mat[RowIndex+i][ColumnIndex+j] = SUB->mat[i][j];
686 }
687 }
688
689 return MAT_OK;
690
691 }
692