← Back to branch summary

~ubuntu-branches/ubuntu/hardy/libterralib/hardy

« back to all changes in this revision

Viewing changes to src/terralib/kernel/TeMatrix.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T Chen
  • Date: 2005-11-25 22:32:59 UTC
  • Revision ID: james.westby@ubuntu.com-20051125223259-3zubal8ux4ki4fjg
Tags: upstream-3.0.3b2
Import upstream version 3.0.3b2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/terralib/kernel/TeMatrix.cpp
 
1
/************************************************************************************
 
2
TerraLib - a library for developing GIS applications.
 
3
Copyright � 2001-2004 INPE and Tecgraf/PUC-Rio.
 
4
 
 
5
This code is part of the TerraLib library.
 
6
This library is free software; you can redistribute it and/or
 
7
modify it under the terms of the GNU Lesser General Public
 
8
License as published by the Free Software Foundation; either
 
9
version 2.1 of the License, or (at your option) any later version.
 
10
 
 
11
You should have received a copy of the GNU Lesser General Public
 
12
License along with this library.
 
13
 
 
14
The authors reassure the license terms regarding the warranties.
 
15
They specifically disclaim any warranties, including, but not limited to,
 
16
the implied warranties of merchantability and fitness for a particular purpose.
 
17
The library provided hereunder is on an "as is" basis, and the authors have no
 
18
obligation to provide maintenance, support, updates, enhancements, or modifications.
 
19
In no event shall INPE and Tecgraf / PUC-Rio be held liable to any party for direct,
 
20
indirect, special, incidental, or consequential damages arising out of the use
 
21
of this library and its documentation.
 
22
*************************************************************************************/
 
23
 
 
24
 // TeMatrix.c -- matrix of type double
 
25
 
 
26
#include <TeMatrix.h>
 
27
#include <math.h>
 
28
#include <stdio.h>
 
29
 
 
30
 
 
31
TeMatrix::TeMatrix(){
 
32
        nrow = 0;
 
33
        ncol = 0;
 
34
        mat  = NULL;
 
35
        lixo = 0.;
 
36
}
 
37
 
 
38
void
 
39
TeMatrix::Clear()
 
40
{
 
41
        if( mat != NULL ){
 
42
                int i;
 
43
                for( i = 0; i < nrow; i++ )
 
44
                        if( mat[i] != NULL )
 
45
                        {
 
46
                                delete []mat[i];
 
47
                                mat[i]=NULL;
 
48
                        }
 
49
                delete []mat;
 
50
                mat = NULL;
 
51
        }
 
52
        nrow = ncol = 0;
 
53
        lixo = (double)0;
 
54
}
 
55
 
 
56
TeMatrix::~TeMatrix()
 
57
{
 
58
        Clear();
 
59
}
 
60
 
 
61
short
 
62
TeMatrix::Alloc(int nl, int nc)
 
63
{
 
64
        Clear();
 
65
        if( nl <= 0 || nc <= 0 ){
 
66
                return false;
 
67
        }
 
68
        nrow = nl;
 
69
        ncol = nc;
 
70
        if (  (mat = new double*[nrow]) == NULL ){
 
71
                Clear();
 
72
                return false;
 
73
        }
 
74
        int i,j;
 
75
        for( i = 0; i < nrow; i++ ) mat[i] = NULL;
 
76
        for( i = 0; i < nrow; i++ ){
 
77
                if( (mat[i] = new double[ncol]) == NULL ){
 
78
                        Clear();
 
79
                        return false;
 
80
                }
 
81
                for( j = 0; j < ncol; j++ ) mat[i][j] = (double)0.;
 
82
        }
 
83
        return true;
 
84
}
 
85
 
 
86
TeMatrix::TeMatrix(const TeMatrix& m)
 
87
{
 
88
     nrow = ncol = 0;
 
89
     mat = NULL;
 
90
     lixo = 0.;
 
91
    *this = m;
 
92
}
 
93
 
 
94
int
 
95
TeMatrix::Init( int nl, int nc , double* f )
 
96
{
 
97
        if( f == NULL )
 
98
                return Init( nl, nc,  (double) 0 );
 
99
        if( Alloc( nl, nc ) ){
 
100
                int c=0;
 
101
                int i,j;
 
102
                for ( i = 0; i < nrow; i++)
 
103
                        for (j = 0; j < ncol; j++)
 
104
                                mat[i][j] = f[c++];
 
105
                return true;
 
106
        }
 
107
        return false;
 
108
}
 
109
 
 
110
int
 
111
TeMatrix::Init(int nl, int nc, double f )
 
112
{
 
113
        // initialize matrix
 
114
        if( Alloc( nl, nc ) ){
 
115
                int i,j;
 
116
                for (i = 0; i < nrow; i++)
 
117
                         for (j = 0; j < ncol; j++)
 
118
                                mat[i][j] = f;
 
119
                return true;
 
120
        }
 
121
        return false;
 
122
}
 
123
 
 
124
int
 
125
TeMatrix::Init(int k, double* f)
 
126
{
 
127
        if( f == NULL )
 
128
                return Init( k, k,  (double) 0 );
 
129
        if( Alloc(k,k) == false ) 
 
130
                return false;
 
131
        int i,j;
 
132
        for ( i = 0; i < nrow; i++  )
 
133
                for (j = 0; j < ncol; j++)
 
134
                        mat[i][j] = (i==j)? f[i]: (double)0.;
 
135
        return true;
 
136
}
 
137
 
 
138
int
 
139
TeMatrix::Init(int k, double f)
 
140
{
 
141
        if( Alloc(k,k) == false )
 
142
                return false;
 
143
        int i,j;
 
144
        for (i=0; i<nrow; i++)
 
145
                for (j=0; j<ncol; j++)
 
146
                        mat[i][j] = (i==j)? (double)f:0;
 
147
        return true;
 
148
}
 
149
 
 
150
TeMatrix& 
 
151
TeMatrix::operator=(const TeMatrix& m)
 
152
{
 
153
 
 
154
        if( Alloc(m.nrow,m.ncol) == false ){
 
155
                printf("\nMatrix::operator=:no memory available \n");
 
156
                return *this;
 
157
        }
 
158
        int     i,j;
 
159
        for( i = 0; i < nrow; i++)
 
160
                for( j = 0; j < ncol; j++)
 
161
                        mat[i][j] = m.mat[i][j];
 
162
        return *this;
 
163
}
 
164
 
 
165
int
 
166
TeMatrix::operator==(const TeMatrix& m) const
 
167
{
 
168
        if( nrow != m.nrow || ncol != m.ncol )
 
169
                return false;
 
170
        int     i,j;
 
171
        for ( i = 0; i < nrow; i++)
 
172
                for ( j = 0; j < ncol; j++)
 
173
                        if ( mat[i][j]!= m.mat[i][j] )
 
174
                                return false;
 
175
        return true;
 
176
}
 
177
 
 
178
void
 
179
TeMatrix::operator*=(double f)
 
180
{
 
181
        int     i,j;
 
182
        for( i = 0; i < nrow; i++)
 
183
                for( j = 0; j < ncol; j++)
 
184
                        mat[i][j] *= f;
 
185
        return;
 
186
}
 
187
 
 
188
 
 
189
//-------------------------------------------------------------friend methods
 
190
 
 
191
TeMatrix
 
192
operator+(const TeMatrix& m,const TeMatrix& n)
 
193
{
 
194
        if( m.nrow != n.nrow || m.ncol != n.ncol ){
 
195
                printf("\nMatrix::operator+ >>> Operandos com tamanhos diferentes\n");
 
196
                return  m;
 
197
        } 
 
198
        TeMatrix rm;
 
199
        if( rm.Init(m.Nrow(),m.Ncol()) == false ){
 
200
                printf("\nMatrix::operator+ >>> Memoria nao disponivel\n");
 
201
                return  m;
 
202
        }
 
203
        int     i,j;
 
204
        for ( i = 0; i < m.Nrow(); i++)
 
205
                for ( j = 0; j < m.Ncol(); j++)
 
206
                        rm.mat[i][j] = m.mat[i][j] + n.mat[i][j];
 
207
        return rm;
 
208
}
 
209
 
 
210
TeMatrix
 
211
operator-(const TeMatrix& m,const TeMatrix& n)
 
212
{
 
213
        if( m.nrow != n.nrow || m.ncol != n.ncol ){
 
214
                printf("\nMatrix::operator+ >>> Operandos com tamanhos diferentes\n");
 
215
                return  m;
 
216
        } 
 
217
        TeMatrix rm;
 
218
        if( rm.Init(m.Nrow(),m.Ncol()) == false ){
 
219
                printf("\nMatrix::operator- >>> Memoria nao disponivel\n");
 
220
                return  m;
 
221
        }
 
222
 
 
223
        int     i,j;
 
224
        for ( i = 0; i <  m.Nrow(); i++)
 
225
                for ( j = 0; j < m.Ncol(); j++)
 
226
                        rm.mat[i][j] = m.mat[i][j] - n.mat[i][j];
 
227
        return rm;
 
228
}
 
229
 
 
230
TeMatrix
 
231
TeMatrix::operator-()
 
232
{
 
233
        TeMatrix rm;  
 
234
        if( rm.Init(Nrow(),Ncol()) == false ){
 
235
                printf("\n operator-:no memory \n");
 
236
                return  *this;
 
237
        }
 
238
        int i,j;
 
239
        for (i=0; i<Nrow(); i++)
 
240
                for (j=0; j<Ncol(); j++)
 
241
                        rm.mat[i][j] = -mat[i][j];
 
242
        return rm;
 
243
}
 
244
 
 
245
TeMatrix 
 
246
operator*(const TeMatrix& m,const TeMatrix& n)
 
247
{
 
248
        if ( m.Ncol() != n.Nrow() ) {
 
249
                printf( "\nMatrix::operator* >>> Operandos com tamanhos diferentes\n");
 
250
                return m;
 
251
        }
 
252
        int nr = m.Nrow(), nc =n.Ncol();
 
253
        TeMatrix result;
 
254
        if( result.Init(nr,nc) == false){
 
255
                printf("\nMatrix::operator* >>> Memoria nao disponivel\n");
 
256
                return  m;
 
257
        }
 
258
 
 
259
        double  sum = (double)0.;
 
260
        int     l,c,i;
 
261
 
 
262
        for ( l = 0; l < m.Nrow(); l++){
 
263
                for ( c = 0; c < n.Ncol(); c++){
 
264
                        sum = (double)0.;
 
265
                        for ( i = 0; i < m.Ncol(); i++)
 
266
                                sum += m.mat[l][i] * n.mat[i][c];
 
267
                        result.mat[l][c] = sum;
 
268
 
 
269
                }
 
270
        }
 
271
        return result; 
 
272
}
 
273
 
 
274
TeMatrix
 
275
operator*(double f,const TeMatrix& m)
 
276
{
 
277
        int nr = m.Nrow(), nc =m.Ncol();
 
278
        TeMatrix rm;
 
279
        if( rm.Init(nr,nc) == false){
 
280
                printf("\noperator*:no memory");
 
281
                return m;
 
282
        }
 
283
        int i,j;
 
284
        for (i = 0; i < rm.Nrow(); i++)
 
285
                for (j = 0; j < rm.Ncol(); j++)
 
286
                        rm.mat[i][j] = f * m.mat[i][j];
 
287
        return rm;
 
288
}
 
289
 
 
290
//VRMC 12/98
 
291
int
 
292
TeMatrix::isUpperTriangle() const
 
293
{
 
294
        // elements above diagonal != 0
 
295
        int i,j;
 
296
        for(j = 0; j < ncol; j++)
 
297
                for(i = j+1; i < nrow; i++)
 
298
                        if ( mat[i][j]  != (double)0 )
 
299
                                return false;
 
300
        return true;
 
301
}
 
302
 
 
303
//VRMC 12/98
 
304
int
 
305
TeMatrix::isLowerTriangle() const       
 
306
{
 
307
        // elements under diagonal different 0
 
308
        int i,j;
 
309
/*      for(j=0; j<ncol; j++)
 
310
                for(i=j+1; i<nrow; i++)
 
311
                        if ( mat[i][j]!=0 )*/
 
312
        // VRMC 11/98
 
313
        for (i=0; i< nrow; i++)
 
314
                for (j= i+1; j <ncol; j++)
 
315
                        if ( mat[i][j] != (double)0)
 
316
                                return false;
 
317
        return true;
 
318
}
 
319
 
 
320
double
 
321
TeMatrix:: Determinant()
 
322
{
 
323
        if (Nrow() != Ncol()) {
 
324
                printf("\nMatrix::Determinant >>> not a square matrix\n");
 
325
                return lixo;
 
326
        }
 
327
        if (Nrow()==1) return mat[0][0];
 
328
        if (Nrow()==2)
 
329
                return mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
 
330
        if (Nrow()==3) 
 
331
                return (  mat[0][0] * mat[1][1] * mat[2][2]
 
332
                         +mat[0][1] * mat[1][2] * mat[2][0]
 
333
                         +mat[0][2] * mat[1][0] * mat[2][1]
 
334
                         -mat[0][1] * mat[1][0] * mat[2][2]
 
335
                         -mat[0][0] * mat[1][2] * mat[2][1]
 
336
                         -mat[0][2] * mat[1][1] * mat[2][0] );
 
337
        if ( isUpperTriangle() ) {
 
338
                double val = 1;
 
339
                int i;
 
340
                for (i = 0; i < Nrow(); i++) val *= mat[i][i];
 
341
                return val;
 
342
        }
 
343
        TeMatrix mt;
 
344
        double val = 0;
 
345
        double det;
 
346
        int sign = 1;
 
347
        int j;
 
348
        for(j = 0; j < Ncol(); j++) {
 
349
                if( CoFactor(0,j, mt) == false )
 
350
                        return (double)0;
 
351
                det = mt.Determinant();         
 
352
                val += sign * mat[0][j] * det;
 
353
                sign *= -1;
 
354
        }
 
355
        return val;
 
356
}
 
357
 
 
358
int
 
359
TeMatrix::CoFactor(int irow, int jcol, TeMatrix& mt ) const
 
360
{
 
361
        if ( nrow == 1||ncol == 1 ) {
 
362
                printf("\nMatrix::CoFactor >>>  can't CoFactor row or column matrix\n");
 
363
                return false;
 
364
        }
 
365
        if( mt.Init(nrow-1,ncol-1) == false){
 
366
                printf("\nMatrix::CoFactor: Memoria nao disponivel\n");
 
367
                return false;
 
368
        }
 
369
        int getcol, getrow =0;
 
370
        int i,j;
 
371
        for(i=0; i < mt.Nrow(); i++) {
 
372
                if ( getrow == irow ) ++getrow;
 
373
                if ( getrow == nrow ) break;
 
374
                getcol = 0;
 
375
                for(j=0; j < mt.Ncol(); j++) {
 
376
                        if ( getcol==jcol ) ++getcol;
 
377
                        if ( getcol==ncol ) continue;
 
378
                        mt.mat[i][j] = mat[getrow][getcol];
 
379
                        ++getcol;
 
380
                }
 
381
                ++getrow;
 
382
        }
 
383
        return true;
 
384
}
 
385
 
 
386
int 
 
387
TeMatrix::Transpose( TeMatrix &rm ) const
 
388
{
 
389
        if( &rm == this ){
 
390
                printf("\nMatrix::Transpose >>> Operacao usa duas matrizes\n");
 
391
                return false;
 
392
        }
 
393
 
 
394
        if( rm.Init( ncol, nrow ) == false){
 
395
                printf("\nMatrix::Transpose >>> Memoria nao disponivel\n");
 
396
                return false;
 
397
        }
 
398
        int i,j;
 
399
        for(i = 0; i < ncol; i++)
 
400
                for( j = 0; j < nrow; j++)
 
401
                        rm.mat[i][j] = mat[j][i];
 
402
        return true;
 
403
}
 
404
 
 
405
// VRMC 12/98
 
406
int
 
407
TeMatrix::isSimetric()  const
 
408
{
 
409
        int i,j;
 
410
 
 
411
        if (nrow != ncol){
 
412
                printf("\nMatrix::isSimetric >>> Memoria nao disponivel\n");
 
413
                return false;
 
414
        }
 
415
 
 
416
        for (i = 0; i < nrow; i++)
 
417
                for (j = 0; j < ncol; j++)
 
418
                        if ( mat[i][j] != mat[j][i]) return false;
 
419
 
 
420
        return  true;
 
421
}
 
422
 
 
423
// VRMC 12/98
 
424
int             
 
425
TeMatrix::isPositiveDefinide() const
 
426
{
 
427
        int i, j, dim, subdim;
 
428
        TeMatrix        Submat;
 
429
 
 
430
        if (nrow != ncol){
 
431
                printf("\nMatrix::isPositiveDefinide >>> Matriz tem que ser quadrada\n");
 
432
                return false;
 
433
        }
 
434
        dim = nrow;
 
435
 
 
436
        for (subdim=1 ; subdim <= dim; subdim++)
 
437
        {
 
438
                if( Submat.Init( subdim,subdim ) == false ){
 
439
                        printf("\nMatrix::isPositiveDefinide>>>Memoria nao disponivel!\n");
 
440
                        return false;
 
441
                }
 
442
 
 
443
                for ( i=0; i < subdim; i++)
 
444
                        for (j = 0; j < subdim; j++)
 
445
                                Submat(i,j) = mat[i][j];
 
446
                
 
447
                if (Submat.Determinant() <= 0.) {
 
448
                        printf("\nMatrix::isPositiveDefinide>>>Matriz nao e positiva definida!\n");
 
449
                        return false;
 
450
                }
 
451
 
 
452
                Submat.Clear();
 
453
        }
 
454
        return true;
 
455
}
 
456
 
 
457
// VRMC 12/98
 
458
int
 
459
TeMatrix :: MatTransf( TeMatrix& mt )
 
460
{
 
461
        TeMatrix        TI;     // inverse of a inferior triangular matrix
 
462
        int     dim;
 
463
 
 
464
        if( &mt == this ){
 
465
                printf("\nMatrix::MatTransf >>> Operacao usa duas matrizes\n");
 
466
                return false;
 
467
        }
 
468
 
 
469
        if( nrow <= 0 || ncol <= 0 ){
 
470
                printf("\nMatrix::MatTransf >>> Dimensoes da matriz invalidas!\n");
 
471
                return false;
 
472
        }
 
473
 
 
474
        if( nrow != ncol ){
 
475
                printf("\nMatrix::MatTransf >>> Dimensoes da matriz invalidas!\n");
 
476
                return false;
 
477
        }
 
478
 
 
479
        dim = nrow;
 
480
        if( mat[0][0] <= 0. ){
 
481
                printf("\nMatrix::MatTranf >>> ERROR\n");
 
482
                return false;
 
483
        }
 
484
 
 
485
        if( TI.Init( dim, (double)0. ) == false ){
 
486
                printf("\nMatrix::MatTransf >>> Memoria nao disponivel!\n");
 
487
                return false;
 
488
        }
 
489
 
 
490
 
 
491
        //--    Calculate inferior triangular matrix by Cholesky decomposition method
 
492
        if ((*this).CholeskyDecomp(TI) == false) {
 
493
                printf("\nMatrix::MatTransf>>> Nao foi possivel decompor a matriz!\n");
 
494
                return false;
 
495
        }
 
496
 
 
497
        //--    Calculate the inverse of Cholesky matrix
 
498
//      if (TI.CholeskyInv(mt) != true)
 
499
        if (!(TI.CholeskyInv(mt)))
 
500
        {
 
501
                printf("\nMatrix::MatTransf>>> Matriz nao inversivel!\n");
 
502
                return false;
 
503
        }
 
504
 
 
505
        return true;
 
506
}
 
507
 
 
508
// VRMC 12/98                   
 
509
int
 
510
TeMatrix::CholeskyInv (TeMatrix&        mt) const
 
511
{
 
512
        if( mt.Init(nrow,(double)0.) == false ){
 
513
                printf("\nMatrix::CholeskyInv >>> Memoria nao disponivel!\n");
 
514
                return false;
 
515
        }
 
516
 
 
517
        if( &mt == this ){
 
518
                printf("\nMatrix::CholeskyInv >>> Operacao usa duas matrizes\n");
 
519
                return false;
 
520
        }
 
521
        if ( nrow != ncol ){
 
522
                printf("\nMatrix::CholeskyInv>>> Can't invert a non-square matrix");
 
523
                return false;
 
524
        }
 
525
//      if (this->isLowerTriangle() != true) {
 
526
        if (!(this->isLowerTriangle())) {
 
527
                printf("\nMatrix::CholeskyInv >>> Matriz nao e triangular inferior\n");
 
528
                return false;
 
529
        }
 
530
 
 
531
        int i; double rf;
 
532
        for( i = 0; i < nrow; i++ ){
 
533
                if( this->mat[i][i] == (double)0. ){
 
534
                        printf("\nMatrix::CholeskyInv >>> ERROR IV\n");
 
535
                        return false;
 
536
                }
 
537
                mt.mat[i][i] = (double)1. / this->mat[i][i];
 
538
                int j,k;
 
539
                for( j = i-1; j >= 0; j-- ){
 
540
                        rf = (double)0.;
 
541
                        for( k = j; k <= i-1; k++ ){
 
542
                                rf += ( this->mat[i][k] * mt.mat[k][j] );
 
543
                        }
 
544
                        mt.mat[i][j] = -rf * mt.mat[i][i];
 
545
                }
 
546
        }
 
547
        return true;
 
548
}
 
549
        
 
550
 
 
551
int
 
552
TeMatrix::EigenValues( TeMatrix& mt )
 
553
{
 
554
        //---   calcula os eigenvalues  de uma
 
555
        //      matriz representada na forma vetorial 
 
556
        //      Os auto valores sao devolvidos em ordem decrescente    
 
557
        //      e ocupam as posicoes:  0, 2, 5, 9, ... .       
 
558
 
 
559
        if( &mt == this ){
 
560
                printf("\nMatrix::EigenValues >>> Operacao usa duas matrizes\n");
 
561
                return false;
 
562
        }
 
563
 
 
564
double  *cov    = NULL,         /* matriz de covariancia        */
 
565
        *e_val  = NULL;         /*   "    de auto_valores       */
 
566
 
 
567
int     dim = (*this).Nrow();           /* ordem da matriz*/
 
568
 
 
569
int     i,j,k,          /*   cov[], e_val[]    */
 
570
        ia,             /*                     */
 
571
        ind,            /*     |0      |       */
 
572
        l,m,            /*     |1 2    |       */
 
573
        mq,lm,          /*     |3 4 5  |       */
 
574
        ll,mm,          /*     |6 7 8 9|       */
 
575
        ilq,imq,        /*                     */
 
576
        im,iq,il,lq,    /*    (para DIM=4)     */
 
577
        jq;
 
578
 
 
579
double  range,          /*       e_vec[]       */
 
580
        dim1,fdim,      /*                     */  
 
581
        anorm,          /*    |0  4  8 12|     */
 
582
        anrmx,          /*    |1  5  9 13|     */
 
583
        thr,            /*    |2  6 10 14|     */
 
584
        x,y,z,          /*    |3  7 11 15|     */
 
585
        sinx,cosx,      /*                     */
 
586
        sinx2,cosx2,
 
587
        sincs;
 
588
 
 
589
int     fat = (dim*dim+dim)/2;
 
590
 
 
591
        if( dim <= 0 ){
 
592
                printf("\nMatrix::EigenValues >>> dimensao da matriz nula!\n");
 
593
                return false;
 
594
        }
 
595
        if( mt.Init( dim, (double)1 ) == false ){
 
596
                printf("\nMatrix::EigenValues >>> Memoria nao disponivel\n");
 
597
                return false;
 
598
        }
 
599
        range = 1.0e-6 ;
 
600
        dim1 = dim - 1 ;
 
601
        fdim = dim;
 
602
 
 
603
        cov = new double[fat];
 
604
 
 
605
        if( cov == NULL ){
 
606
                printf("\nMatrix::EigenValues >>> Memoria nao disponivel (COV)\n");
 
607
                return false;
 
608
        }
 
609
 
 
610
        e_val = new double[fat];
 
611
        if( e_val == NULL ){
 
612
                delete []cov; //SSL0296
 
613
                cov = NULL;
 
614
                printf("\nMatrix::EigenValues >>> Memoria nao disponivel (EVAL)\n");
 
615
                return false;
 
616
        }
 
617
        k = 0;
 
618
        for ( i = 0; i< dim; i++)
 
619
                for ( j = 0; j <= i; j++)
 
620
                        cov[k++] = mat[i][j];
 
621
 
 
622
        for (i=0; i < ((dim*dim+dim)/2); i++) e_val[i] = cov[i];
 
623
 
 
624
        anorm = 0.0;
 
625
 
 
626
        for ( j = 0; j < dim; j++){
 
627
                for (i = 0; i <= j; i++){
 
628
                        if (i != j){
 
629
                                ia = i + (j*j+j)/2;
 
630
                                anorm = anorm + e_val[ia] * e_val[ia];
 
631
                        }
 
632
                }
 
633
        }
 
634
 
 
635
        if (anorm <= 0) goto l5;
 
636
 
 
637
        anorm = 1.414 * sqrt((double)anorm);
 
638
        anrmx = anorm * range/fdim;
 
639
        ind   = 0;
 
640
        thr   = anorm;
 
641
 
 
642
l1:
 
643
        thr = thr/fdim;
 
644
l2:
 
645
        l = 0;
 
646
l3:
 
647
        m = l+1;
 
648
l4:
 
649
        mq = (m*m + m)/2;
 
650
        lq = (l*l + l)/2;
 
651
        lm = l + mq;
 
652
 
 
653
        if ( fabs((double)(e_val[lm])) >= thr )    
 
654
        {
 
655
                ind = 1;
 
656
                ll = l + lq;
 
657
                mm = m + mq;
 
658
 
 
659
                x = 0.5 * (e_val[ll] - e_val[mm]);
 
660
                z = e_val[lm] * e_val[lm] + x*x;
 
661
                y = - e_val[lm] / sqrt((double)(z));
 
662
 
 
663
                if (x < 0)  
 
664
                {
 
665
                        y = -y;
 
666
                }
 
667
 
 
668
                z = sqrt( (double)(1.0-y*y) );
 
669
                sinx = y / sqrt( (double)( 2.0*(1.0 + z) ) );
 
670
                sinx2 = sinx * sinx;
 
671
                
 
672
                cosx = sqrt( (double)(1.0 - sinx2) );
 
673
                cosx2 = cosx * cosx;
 
674
 
 
675
                sincs = sinx * cosx;
 
676
        
 
677
                ilq = dim * l;
 
678
                imq = dim * m;
 
679
 
 
680
                for (i = 0; i < dim; i++)
 
681
                {
 
682
                        iq = (i*i + i)/2;
 
683
                        if ( i != l )
 
684
                        {
 
685
                                if (i != m)
 
686
                                {
 
687
                                        if (i > m) 
 
688
                                                im = m + iq;
 
689
                                        else
 
690
                                                im = i + mq;
 
691
                                        if (i < l)
 
692
                                                il = i + lq;
 
693
                                        else
 
694
                                                il = l + iq;
 
695
 
 
696
                                        x = e_val[il] * cosx - e_val[im] * sinx;
 
697
                                        e_val[im] = e_val[il] * sinx + e_val[im] * cosx;
 
698
                                        e_val[il] = x;
 
699
                                }
 
700
                        }
 
701
                }
 
702
                x = 2.0 * e_val[lm] * sincs;
 
703
                y = e_val[ll] * cosx2 + e_val[mm] * sinx2 - x;
 
704
                x = e_val[ll] * sinx2 + e_val[mm] * cosx2 + x;
 
705
 
 
706
                e_val[lm] = (e_val[ll]-e_val[mm])*sincs+e_val[lm]*(cosx2-sinx2);
 
707
                e_val[ll] = y;
 
708
                e_val[mm] = x;
 
709
        }
 
710
 
 
711
        if (m != (dim-1))
 
712
        {
 
713
                m = m + 1;      goto l4;
 
714
        }
 
715
        if (l != (dim-2))
 
716
        {
 
717
                l = l + 1;      goto l3;
 
718
        }
 
719
        if (ind == 1)
 
720
        {
 
721
                ind = 0;        goto l2;
 
722
        }
 
723
        if (thr > anrmx) goto l1;
 
724
 
 
725
l5:
 
726
        iq = -dim;
 
727
 
 
728
        for (i = 0; i < dim; i++)
 
729
        {
 
730
                iq = iq + dim;
 
731
                ll = i + (i*i + i)/2;
 
732
                jq = dim * (i-1);
 
733
 
 
734
                for (j = i; j < dim; j++)
 
735
                {
 
736
                        jq = jq + dim;
 
737
                        mm = j + (j*j + j)/2;
 
738
 
 
739
                        if (e_val[ll] < e_val[mm])
 
740
                        {
 
741
                                x = e_val[ll];
 
742
                                e_val[ll] = e_val[mm];
 
743
                                e_val[mm] = x;
 
744
                        }
 
745
                }
 
746
        }
 
747
 
 
748
        for ( i = 0; i< dim; i++){
 
749
                mt(i,0) = e_val[(i*(i+1))/2+i];
 
750
        }
 
751
 
 
752
        delete []cov;   //SSL0296
 
753
        delete []e_val; //SSL0296
 
754
 
 
755
        return true;
 
756
}
 
757
 
 
758
//--------------------------------------------------EigenVectors
 
759
//
 
760
//      Metodo adaptado a partir de uma implementacao em C
 
761
//      utilizada no software SITIM.
 
762
 
 
763
int
 
764
TeMatrix::EigenVectors( TeMatrix& mt )
 
765
{
 
766
        if( &mt == this ){
 
767
                printf("\nMatrix::EigenVectors >>> Operacao usa duas matrizes\n");
 
768
                return false;
 
769
        }
 
770
 
 
771
double  *cov    = NULL,         /* matriz de covariancia        */
 
772
        *e_val  = NULL,         /*   "    de auto_valores       */
 
773
        *e_vec  = NULL;         /*   "    de auto_vetores       */
 
774
 
 
775
int     dim = (*this).Nrow();           /* ordem da matriz*/
 
776
 
 
777
int     i,j,ij,k,       /*   cov[], e_val[]    */
 
778
        ia,             /*                     */
 
779
        ind,            /*     |0      |       */
 
780
        l,m,            /*     |1 2    |       */
 
781
        mq,lm,          /*     |3 4 5  |       */
 
782
        ll,mm,          /*     |6 7 8 9|       */
 
783
        ilq,imq,        /*                     */
 
784
        im,iq,il,lq,    /*    (para DIM=4)     */
 
785
        ilr,imr,
 
786
        jq;
 
787
 
 
788
double  range,          /*       e_vec[]       */
 
789
        dim1,fdim,      /*                     */  
 
790
        anorm,          /*    |0  4  8 12|     */
 
791
        anrmx,          /*    |1  5  9 13|     */
 
792
        thr,            /*    |2  6 10 14|     */
 
793
        x,y,z,          /*    |3  7 11 15|     */
 
794
        sinx,cosx,      /*                     */
 
795
        sinx2,cosx2,
 
796
        sincs;
 
797
 
 
798
        if( dim <= 0 ){
 
799
                printf("\nMatrix::EigenValues >>> dimensao da matriz nula!\n");
 
800
                return false;
 
801
        }
 
802
        if( mt.Init( dim, (double)1 ) == false ){
 
803
                printf("\nMatriz::EigenVectors >>> Memoria nao disponivel!\n");
 
804
                return false;
 
805
        }
 
806
 
 
807
        int fat =(dim*dim+dim)/2;
 
808
        range = 1.0e-6 ;
 
809
        dim1 = dim - 1 ;
 
810
        fdim = dim;
 
811
 
 
812
        cov   = new double[fat];
 
813
        e_vec = new double[dim*dim];
 
814
        e_val = new double[fat];
 
815
 
 
816
        if( cov == NULL || e_vec == NULL || e_val == NULL ){
 
817
                printf("\nMatrix::EigenVectors >>> Memoria nao disponivel\n");
 
818
                return false;
 
819
        }
 
820
 
 
821
        k = 0;
 
822
        for ( i = 0; i< dim; i++)
 
823
                for ( j = 0; j <= i; j++)
 
824
                        cov[k++] = (*this)(i,j);
 
825
 
 
826
        for (i=0; i < ((dim*dim+dim)/2); i++) e_val[i] = cov[i];
 
827
 
 
828
        iq = -dim;
 
829
        for (i = 0; i < dim; i++)
 
830
        {
 
831
                iq = iq + dim;
 
832
                for ( j = 0; j < dim; j++)
 
833
                {
 
834
                        ij = iq + j;
 
835
                        if (i == j)
 
836
                                e_vec[ij] = 1.0;
 
837
                        else
 
838
                                e_vec[ij] = 0.0;
 
839
                }
 
840
        }
 
841
 
 
842
        anorm = 0.0;
 
843
 
 
844
        for ( j = 0; j < dim; j++)
 
845
        {       
 
846
                for (i = 0; i <= j; i++)
 
847
                {
 
848
                        if (i != j)
 
849
                        {
 
850
                                ia = i + (j*j+j)/2;
 
851
                                anorm = anorm + e_val[ia] * e_val[ia];
 
852
                        }
 
853
                }
 
854
        }
 
855
 
 
856
        if (anorm <= 0) goto l5;
 
857
 
 
858
        anorm = 1.414 * sqrt((double)anorm);
 
859
        anrmx = anorm * range/fdim;
 
860
        ind   = 0;
 
861
        thr   = anorm;
 
862
 
 
863
l1:
 
864
        thr = thr/fdim;
 
865
l2:
 
866
        l = 0;
 
867
l3:
 
868
        m = l+1;
 
869
l4:
 
870
        mq = (m*m + m)/2;
 
871
        lq = (l*l + l)/2;
 
872
        lm = l + mq;
 
873
 
 
874
        if ( fabs((double)(e_val[lm])) >= thr )    
 
875
        {
 
876
                ind = 1;
 
877
                ll = l + lq;
 
878
                mm = m + mq;
 
879
 
 
880
                x = 0.5 * (e_val[ll] - e_val[mm]);
 
881
                z = e_val[lm] * e_val[lm] + x*x;
 
882
                y = - e_val[lm] / sqrt((double)(z));
 
883
 
 
884
                if (x < 0)  
 
885
                {
 
886
                        y = -y;
 
887
                }
 
888
 
 
889
                z = sqrt( (double)(1.0-y*y) );
 
890
                sinx = y / sqrt( (double)( 2.0*(1.0 + z) ) );
 
891
                sinx2 = sinx * sinx;
 
892
                
 
893
                cosx = sqrt( (double)(1.0 - sinx2) );
 
894
                cosx2 = cosx * cosx;
 
895
 
 
896
                sincs = sinx * cosx;
 
897
        
 
898
                ilq = dim * l;
 
899
                imq = dim * m;
 
900
 
 
901
                for (i = 0; i < dim; i++)
 
902
                {
 
903
                        iq = (i*i + i)/2;
 
904
                        if ( i != l )
 
905
                        {
 
906
                                if (i != m)
 
907
                                {
 
908
                                        if (i > m) 
 
909
                                                im = m + iq;
 
910
                                        else
 
911
                                                im = i + mq;
 
912
                                        if (i < l)
 
913
                                                il = i + lq;
 
914
                                        else
 
915
                                                il = l + iq;
 
916
 
 
917
                                        x = e_val[il] * cosx - e_val[im] * sinx;
 
918
                                        e_val[im] = e_val[il] * sinx + e_val[im] * cosx;
 
919
                                        e_val[il] = x;
 
920
                                }
 
921
                        }
 
922
 
 
923
                        //---   calculate eigenvectors
 
924
 
 
925
                        ilr = ilq + i;
 
926
                        imr = imq + i;
 
927
                        x   = e_vec[ilr] * cosx - e_vec[imr] * sinx;
 
928
                        e_vec[imr] = e_vec[ilr] * sinx + e_vec[imr] * cosx;
 
929
                        e_vec[ilr] = x;
 
930
                }
 
931
 
 
932
                x = 2.0 * e_val[lm] * sincs;
 
933
                y = e_val[ll] * cosx2 + e_val[mm] * sinx2 - x;
 
934
                x = e_val[ll] * sinx2 + e_val[mm] * cosx2 + x;
 
935
 
 
936
                e_val[lm] = (e_val[ll]-e_val[mm])*sincs+e_val[lm]*(cosx2-sinx2);
 
937
                e_val[ll] = y;
 
938
                e_val[mm] = x;
 
939
        }
 
940
 
 
941
        if (m != (dim-1))
 
942
        {
 
943
                m = m + 1;      goto l4;
 
944
        }
 
945
        if (l != (dim-2))
 
946
        {
 
947
                l = l + 1;      goto l3;
 
948
        }
 
949
        if (ind == 1)
 
950
        {
 
951
                ind = 0;        goto l2;
 
952
        }
 
953
        if (thr > anrmx) goto l1;
 
954
 
 
955
l5:
 
956
        iq = -dim;
 
957
 
 
958
        for (i = 0; i < dim; i++)
 
959
        {
 
960
                iq = iq + dim;
 
961
                ll = i + (i*i + i)/2;
 
962
                jq = dim * (i-1);
 
963
 
 
964
                for (j = i; j < dim; j++)
 
965
                {
 
966
                        jq = jq + dim;
 
967
                        mm = j + (j*j + j)/2;
 
968
 
 
969
                        if (e_val[ll] < e_val[mm])
 
970
                        {
 
971
                                x = e_val[ll];
 
972
                                e_val[ll] = e_val[mm];
 
973
                                e_val[mm] = x;
 
974
 
 
975
                                for (k = 0; k < dim; k++)
 
976
                                {
 
977
                                        ilr = iq + k;
 
978
                                        imr = jq + k;
 
979
                                        x   = e_vec[ilr];
 
980
                                        e_vec[ilr] = e_vec[imr];
 
981
                                        e_vec[imr] = x;
 
982
                                }
 
983
                        }
 
984
                }
 
985
        }
 
986
        
 
987
        k=0;
 
988
        for ( i = 0; i< dim; i++){
 
989
                for ( j = 0; j< dim; j++)
 
990
                        mt(j,i) = e_vec[k++];
 
991
        }
 
992
 
 
993
        delete []cov;   //SSL0296
 
994
        delete []e_vec; //SSL0296
 
995
        delete []e_val; //SSL0296
 
996
 
 
997
        return true;
 
998
}
 
999
 
 
1000
//--------------------------------------------------EigenVec
 
1001
//
 
1002
//      Metodo desenvolvido para suporte de decisao (suporte_stubs.cpp)
 
1003
//      Missae Yamamoto  (junho/1999)
 
1004
 
 
1005
int
 
1006
TeMatrix::EigenVec( double e_vec[] )
 
1007
{
 
1008
 
 
1009
double  *e_vec_aux      = NULL, soma = 0.0; 
 
1010
        
 
1011
int     dim = (*this).Nrow();           /* ordem da matriz*/
 
1012
 
 
1013
TeMatrix  aux1, aux2, mt;
 
1014
 
 
1015
int     i,j,k;       
 
1016
        
 
1017
        if( dim <= 0 )
 
1018
        {
 
1019
                printf("\nMatrix::EigenVecVal >>> dimensao da matriz nula!\n");
 
1020
                return false;
 
1021
        }
 
1022
 
 
1023
        if( aux1.Init( dim, (double)1 ) == false )
 
1024
        {
 
1025
                printf("\nMatrix::EigenVecVal >>> Memoria nao disponivel!\n");
 
1026
                return false;
 
1027
        }
 
1028
 
 
1029
        if( aux2.Init( dim, (double)1 ) == false )
 
1030
        {
 
1031
                printf("\nMatrix::EigenVecVal >>> Memoria nao disponivel!\n");
 
1032
                return false;
 
1033
        }
 
1034
 
 
1035
        e_vec_aux = new double[dim];
 
1036
 
 
1037
        if( e_vec_aux == NULL )
 
1038
        {
 
1039
                printf("\nMatrix::EigenVecVal >>> Memoria nao disponivel\n");
 
1040
                return false;
 
1041
        }
 
1042
 
 
1043
        for ( k = 0; k< dim; k++)
 
1044
        {
 
1045
                e_vec_aux[k] = 0.0;
 
1046
                e_vec[k] = 0.0;
 
1047
        }
 
1048
 
 
1049
        aux1 = *this;
 
1050
        aux2 = *this;
 
1051
        
 
1052
        for (;;)
 
1053
        {
 
1054
                mt = aux1 * aux2;
 
1055
 
 
1056
                for (i=0; i<dim; i++)
 
1057
                {
 
1058
                        for (j=0; j<dim; j++)
 
1059
                                 e_vec[i] = e_vec[i] + mt(i,j);
 
1060
                        soma = soma + e_vec[i];
 
1061
                }
 
1062
 
 
1063
                for (j=0; j<dim; j++)
 
1064
                        e_vec[j] = e_vec[j] / soma;
 
1065
 
 
1066
                for (j=0; j<dim; j++)
 
1067
                {
 
1068
                        if ( fabs(e_vec_aux[j] - e_vec[j]) < 0.001 ) 
 
1069
                        {
 
1070
                                delete []e_vec_aux;     
 
1071
                                return true;
 
1072
                        }
 
1073
                        e_vec_aux[j] = e_vec[j];
 
1074
                }
 
1075
 
 
1076
                aux1 = mt;
 
1077
                aux2 = mt;
 
1078
        }
 
1079
 
 
1080
}
 
1081
 
 
1082
//------------------------------------------ Print
 
1083
void
 
1084
TeMatrix::Print()
 
1085
{
 
1086
        //printf("\nMATRIZ ( %d, %d ) LIXO %f **mat %p\n", nrow,ncol, lixo, mat);
 
1087
        printf("\nMATRIZ ( %d, %d )\n\n", nrow,ncol);
 
1088
        if( mat == NULL ){
 
1089
                printf("\n>>> mat e NULL \n");
 
1090
                return;
 
1091
        }
 
1092
        int i,j;        
 
1093
        for(i = 0; i < nrow; i++){
 
1094
                for( j = 0; j < ncol;j++){
 
1095
                        printf("%10.4f ",(float)mat[i][j]);
 
1096
                }
 
1097
                printf("\n\n");
 
1098
        }
 
1099
}
 
1100
 
 
1101
// VRMC 12/98
 
1102
int
 
1103
TeMatrix::CholeskyDecomp( TeMatrix& mt )
 
1104
{
 
1105
 
 
1106
        //--    Verify if the matrix is simetric
 
1107
        if ((*this).isSimetric() == false){
 
1108
                printf("\nMatrix::CholeskyDecomp >>> Matriz nao e simetrica!\n");
 
1109
                return false;
 
1110
        }
 
1111
        //--    Verify if the matrix is positive definide
 
1112
        if ((*this).isPositiveDefinide() == false){
 
1113
                printf("\nMatrix::CholeskyDecomp >>> Matriz nao e positiva definida!\n");
 
1114
                return false;
 
1115
        }
 
1116
        if( &mt == this ){
 
1117
                printf("\nMatrix::CholeskyDecomp >>> Operacao usa duas matrizes\n");
 
1118
                return false;
 
1119
        }
 
1120
        if( mat[0][0] <= 0. ){
 
1121
                printf("\nMatrix::CholeskyDecomp >>> Posicao [0,0] e' nula\n");
 
1122
                return false;
 
1123
        }
 
1124
    if( mt.Init(nrow,ncol) == false ){
 
1125
                printf("\nMatrix::CholeskyDecomp>>>no memory \n");
 
1126
                return false;
 
1127
        }
 
1128
        mt(0,0) = (double)sqrt((double)(*this)(0,0));
 
1129
 
 
1130
        //---   Calculate the first column of the TeMatrix mt
 
1131
 
 
1132
        int i,j,k;
 
1133
        for( i= 0; i < nrow; i++){
 
1134
                mt(i,0) = (*this)(i,0)/mt(0,0);
 
1135
        }
 
1136
         
 
1137
        for( i= 1; i < nrow; i++){
 
1138
                for( j= 1; j <= i; j++){
 
1139
                        //int m = (j*j - j)/2 + 1;
 
1140
                        double rf = 0.;
 
1141
                        //int k2 = j - 1;
 
1142
                        //for( k = 0;k<=k2;k++)
 
1143
                        for( k = 0;k < j;k++)
 
1144
                                rf += mt(i,k) * mt(j,k);
 
1145
                        if( i == j ){
 
1146
                                rf = (*this)(i,j) - rf;
 
1147
                                if( rf < 0. ){
 
1148
                                        printf("\nMatrix::CholeskyDecomp:ERRO \n");
 
1149
                                        return false;
 
1150
                                }
 
1151
                                mt(i,j) = (double)sqrt((double)rf);
 
1152
                        }
 
1153
                        else{
 
1154
                                if( mt(j,j) == 0. ){
 
1155
                                        printf("\nMatrix::CholeskyDecomp:ERRO \n");
 
1156
                                        return false;
 
1157
                                }
 
1158
                                mt(i,j) = ((*this)(i,j)-rf)/mt (j,j );
 
1159
                        }
 
1160
                }
 
1161
        }
 
1162
        return true;
 
1163
}
 
1164
 
 
1165
 
 
1166
//-------------------------------------------------Inverse(SITIM)
 
1167
//
 
1168
//---   Calcula a matriz inversa pelo metodo de Gauss-Jordan.
 
1169
//
 
1170
 
 
1171
int
 
1172
TeMatrix::Inverse( TeMatrix& matinv ) const
 
1173
{
 
1174
 
 
1175
        if( &matinv == this ){
 
1176
                printf("\nMatrix::Inverse >>> Operacao usa duas matrizes\n");
 
1177
                return false;
 
1178
        }
 
1179
        if ( nrow != ncol ){
 
1180
                printf("\nMatrix::Inverse: can't invert a non-square matrix");
 
1181
                return false;
 
1182
         }
 
1183
 
 
1184
        if( matinv.Init( this->nrow, (double)1 ) == false ){
 
1185
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1186
                return false;
 
1187
        }
 
1188
 
 
1189
//---   Calculo da inversa se a matriz tem dimensao unitaria
 
1190
 
 
1191
        if(nrow == 1){
 
1192
                if(mat[0][0] == 0.0) {
 
1193
                        printf("\noperator~: can't invert singular matrix");
 
1194
                        return false;
 
1195
                }
 
1196
                matinv(0,0) = 1. / mat[0][0];
 
1197
                return true;
 
1198
        }
 
1199
  
 
1200
//---   Formacao da matriz aumentada mataum
 
1201
 
 
1202
        int     m2   = 2 * nrow ;
 
1203
        TeMatrix        mataum;
 
1204
 
 
1205
        if( mataum.Init( nrow , m2 ) == false ){
 
1206
           printf("\nInverse:no memory");
 
1207
           return false;
 
1208
        }
 
1209
        
 
1210
        int i1,i2,i3;
 
1211
        for(i1 = 0 ; i1 < nrow ; i1++){
 
1212
                for( i2 = 0 ; i2 < ncol ; i2++)
 
1213
                        mataum(i1,i2) = mat[i1][i2];
 
1214
 
 
1215
                for( i3 = nrow ; i3 < m2 ; i3++)
 
1216
                        mataum(i1,i3) = 0.0;
 
1217
 
 
1218
                mataum( i1, i1 + nrow) = 1.0;
 
1219
        }
 
1220
 
 
1221
//--    Inicializa ponteiro de linhas
 
1222
 
 
1223
        double  *maxlinha = NULL;
 
1224
 
 
1225
        if( ( maxlinha = new double[nrow] ) == NULL ){
 
1226
                printf("\nMatrix::Inverse>>no memory");
 
1227
                return false;
 
1228
        }
 
1229
 
 
1230
        double  sup,lcon1,amul,pivo;
 
1231
        int     * lp = NULL;
 
1232
 
 
1233
        if( (lp = new int[nrow]) == NULL ){
 
1234
                printf("\nMatrix::Inverse>>no memory");
 
1235
                if( maxlinha != NULL){delete []maxlinha; maxlinha = NULL;} //SSL0296
 
1236
                return false;
 
1237
        }
 
1238
 
 
1239
        for(i1 = 0 ; i1 < nrow ; i1++)
 
1240
                lp[i1] = i1;
 
1241
 
 
1242
        int     lcon = 0;
 
1243
        int     m1 = nrow - 1;
 
1244
 
 
1245
        while(lcon < nrow){
 
1246
//---           Selecao do maior valor em cada linha . seu inverso e armazenado em maxlinha.
 
1247
                for(i1 = 0 ; i1 < nrow ; i1++){
 
1248
                        sup = 0.0;
 
1249
                        for(i2 = 0 ; i2 < nrow ; i2++){
 
1250
                                lcon1 = (double)fabs(mataum(lp[i1],i2));
 
1251
                                if((lcon1 - sup ) > 0.000001 )
 
1252
                                sup = lcon1;
 
1253
                        }
 
1254
                        if(sup == 0.0){
 
1255
                                printf("\nMatrix::Inverse>>Determinant = 0.");
 
1256
                                if( maxlinha != NULL){ delete []maxlinha; maxlinha = NULL; } //SSL0296
 
1257
                                if( lp != NULL){ delete []lp; lp = NULL; } //SSL0296
 
1258
                                return false;
 
1259
                        }
 
1260
                        maxlinha[lp[i1]] = 1.0 / sup;
 
1261
                 }
 
1262
 
 
1263
//---           Selecao do pivo em cada linha
 
1264
 
 
1265
                int lpi;
 
1266
                double  supc = fabs((double)( mataum(lp[lcon],lcon) * maxlinha[lp[lcon]] ));
 
1267
                pivo = mataum(lp[lcon],lcon);
 
1268
                lpi  = lcon;
 
1269
                for(i1 = lcon ; i1 < nrow ; i1++){
 
1270
                        amul = mataum(lp[i1],lcon) * maxlinha[lp[i1]];
 
1271
                        if (fabs(amul) > supc ){
 
1272
                        supc = (fabs((double)amul));
 
1273
                        pivo = mataum(lp[i1],lcon);
 
1274
                        lpi  = i1;
 
1275
                        }
 
1276
                }
 
1277
 
 
1278
                if(pivo == 0.0){
 
1279
                        printf("\nMatrix::Inverse>>Determinant = 0.");
 
1280
                        if( maxlinha != NULL){ delete []maxlinha; maxlinha = NULL; } //SSL0296
 
1281
                        if( lp != NULL){ delete []lp; lp = NULL; } //SSL0296
 
1282
                        return false;
 
1283
                }
 
1284
     
 
1285
                /* troca de linhas */
 
1286
                i1       = lp[lcon];
 
1287
                lp[lcon] = lp[lpi];
 
1288
                lp[lpi]  = i1;
 
1289
        
 
1290
//---           Divide a linha pelo pivo
 
1291
                for(i2 = 0 ; i2 < m2 ; i2++){
 
1292
                        mataum(lp[lcon],i2) = mataum(lp[lcon],i2) / pivo ;
 
1293
                }
 
1294
 
 
1295
//---           Realiza subtracoes nas linhas de forma a zerar os elementos da coluna
 
1296
//---           dada pelo indice lcon
 
1297
 
 
1298
                int     il = lcon;
 
1299
 
 
1300
                for(i2 = 0; i2 < m1 ; i2++){
 
1301
                        il = il + 1;
 
1302
                        if(il >= nrow ) il = 0;
 
1303
                        double  aux = - mataum(lp[il],lcon);
 
1304
                        for(i3 = lcon; i3 < m2 ; i3++)
 
1305
                                mataum(lp[il],i3) = mataum(lp[il],i3) + aux * mataum(lp[lcon],i3);
 
1306
 
 
1307
                        mataum(lp[il],lcon) = 0.0;
 
1308
                }
 
1309
 
 
1310
                lcon = lcon + 1;
 
1311
        }
 
1312
 
 
1313
//---   Copia a parte extendida de mataum em matinv que passa a ser a matriz inversa
 
1314
        for(i1 = 0 ; i1 < nrow ; i1++){
 
1315
                for(i2 = nrow ; i2 < m2 ; i2++){
 
1316
                        i3 = i2 - nrow ;
 
1317
                        matinv(i1,i3) = mataum(lp[i1],i2);
 
1318
                }
 
1319
        }
 
1320
 
 
1321
        if( maxlinha != NULL){
 
1322
                delete []maxlinha; //SSL0296
 
1323
                maxlinha = NULL;
 
1324
        }
 
1325
        if( lp != NULL){
 
1326
                delete []lp; //SSL0296
 
1327
                lp = NULL;
 
1328
        }
 
1329
 
 
1330
        return true;
 
1331
}
 
1332
 
 
1333
//      Esse m�todo, n�o esta correto 
 
1334
//------------------------------------------Inverse
 
1335
//      NIH library
 
1336
// 1. triangulate: *this = ~P*T
 
1337
// 2. when T isUpperTriangle ~T isUpperTriangle
 
1338
// 3. split: T = T.row(0) + subT
 
1339
// 4. I = ~T*T
 
1340
// 5. ~T.at(0,0) = 1/T.at(0,0)
 
1341
// 6. sub~T = ~(subT)
 
1342
// 7. ~T.row(0) = [1/T.at(0,0)]&B
 
1343
//    where T.at(0,0)*B = [t21 ... t2n]*~subT
 
1344
// 8. ~*this = ~T*P
 
1345
 
 
1346
/*int 
 
1347
TeMatrix::InverseNIH( TeMatrix& mt ) const
 
1348
{
 
1349
        TeMatrix  r;
 
1350
        TeMatrix        _subT;
 
1351
        TeMatrix        B;
 
1352
        TeMatrix        val;
 
1353
 
 
1354
        if ( nrow!=ncol || nrow <= 0 || ncol <= 0 ){
 
1355
                printf("\nMatrix::Inverse >>> Matriz nao inversivel!\n");
 
1356
                return false;
 
1357
        }
 
1358
 
 
1359
        if( mt.Init( this->nrow, (double)1 ) == false ){
 
1360
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1361
                return false;
 
1362
        }
 
1363
 
 
1364
        if ( nrow==1 ){
 
1365
                if( mat[0][0] == (double)0 ){
 
1366
                        printf("\nMatrix::Inverse >>> Matriz nao inversivel!\n");
 
1367
                        return false;
 
1368
                }
 
1369
                mt.mat[0][0] = 1/mat[0][0];
 
1370
                return true;
 
1371
        }
 
1372
        TeMatrix T(*this);
 
1373
        if( T.Initialized() == false ){
 
1374
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1375
                return false;
 
1376
        }
 
1377
                
 
1378
        // 1. triangulate: *this = ~P*T
 
1379
        TeMatrix P;
 
1380
        if( P.Init(nrow,ncol) == false){
 
1381
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1382
                return false;
 
1383
        }
 
1384
        if( T.UpperTriangle( P ) == false ){
 
1385
                printf("\nMatriz::Inverse >>> ERRO!\n");
 
1386
                return false;
 
1387
        }
 
1388
        if ( T.Determinant()==0 ) {
 
1389
                printf("\nMatrix::Inverse >>> can't invert singular matrix\n");
 
1390
                return false;
 
1391
        }      
 
1392
 
 
1393
        // 2. when T isUpperTriangle ~T isUpperTriangle
 
1394
        // 3. split: T = T.row(0) + subT
 
1395
 
 
1396
        if( r.Init(1,ncol-1,&(T.mat[0][1])) == false){
 
1397
                printf("\nMatrix::Inverse* >>> Memoria nao disponivel\n");
 
1398
                return false;
 
1399
        }
 
1400
 
 
1401
        T.CoFactor(0,0,mt);
 
1402
        mt.Inverse(_subT);
 
1403
        if( _subT.Initialized()  == false){
 
1404
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1405
                return false;
 
1406
        }
 
1407
 
 
1408
 
 
1409
        B = (-(1/T.mat[0][0])*r*_subT);
 
1410
        if( B.Initialized() == false ){
 
1411
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1412
                return false;
 
1413
        }
 
1414
 
 
1415
        if( val.Init(nrow,ncol) == false){
 
1416
                printf("\nMatrix::Inverse >>> Memoria nao disponivel\n");
 
1417
                return false;
 
1418
        }
 
1419
 
 
1420
        val.mat[0][0] = 1/T.mat[0][0];
 
1421
        int i;
 
1422
        for(i=1; i<ncol; i++) val.mat[0][i]=B.mat[0][i-1];
 
1423
        val.CoFactor(0,0,_subT);
 
1424
 
 
1425
        mt = val*P; // P is now the row-reduction transformation
 
1426
 
 
1427
        return true; 
 
1428
}
 
1429
*/
 
1430
 
 
1431
// os m�todos switchRows e combineRows ficarao comentados 
 
1432
// pois s� sao utilizados no calculo da matriz inversa NIH,
 
1433
// inversaNIH, que tambem foi comentado 
 
1434
/*void 
 
1435
TeMatrix::switchRows(int i,int j)
 
1436
{
 
1437
        double *tmp;
 
1438
 
 
1439
        if( i < 0 || j < 0 || i >= nrow || j >= ncol ){
 
1440
                printf("\nMatrix::switchRows >>> Switching linhas invalidas\n");
 
1441
                return;
 
1442
        }
 
1443
        tmp    = mat[i];
 
1444
        mat[i] = mat[j];
 
1445
        mat[j] = tmp; 
 
1446
}
 
1447
 
 
1448
void
 
1449
TeMatrix::combineRows(int i, double a, int j)
 
1450
{
 
1451
        if( i < 0 || j < 0 || i >= nrow || j >= ncol ){
 
1452
                printf("\nMatrix::combineRows >>> Combining linhas invalidas\n");
 
1453
                return;
 
1454
        }
 
1455
        int h;
 
1456
        for(h = 0; h < ncol; h++)
 
1457
        mat[i][h] += a * mat[j][h];
 
1458
}
 
1459
*/
 
1460
 
 
1461
//--------Esse metodo nao esta correto e por isso ficar� comentado VRMC
 
1462
//----------------------------------------------------
 
1463
/*int 
 
1464
TeMatrix::UpperTriangle( TeMatrix& mt )
 
1465
{
 
1466
        if( &mt == this ){
 
1467
                printf("\nMatrix::UpperTriangle >>> Operacao usa duas matrizes\n");
 
1468
                return false;
 
1469
        }
 
1470
        if( mt.Init(nrow,(double)1) == false){
 
1471
                printf("\nMatrix::UpperTriangle >>> Memoria nao disponivel\n");
 
1472
                return  false;
 
1473
        }       
 
1474
 
 
1475
        if ( this->isUpperTriangle() )
 
1476
                return true;
 
1477
        
 
1478
        int j;
 
1479
        for(j = 0; j < ncol; j++) {
 
1480
                int b_row = nrow-1;  // 1st non-zero entry from bottom
 
1481
                int t_row = j;       // 1st zero entry from the top
 
1482
 
 
1483
                //--- switch rows until all zeros are at
 
1484
                //--- the bottom of jTH column 
 
1485
 
 
1486
                while ( b_row>=t_row ) { 
 
1487
                        while ( b_row > j && mat[b_row][j] == 0. )
 
1488
                                --b_row;
 
1489
                        if ( b_row == j )
 
1490
                                break; // bottom at diagonal
 
1491
                        while ( b_row >= t_row && mat[t_row][j] != 0. )
 
1492
                                ++t_row;
 
1493
                        if ( t_row == Nrow() )
 
1494
                                break; // top at last row
 
1495
                        if ( t_row > b_row )
 
1496
                                break;
 
1497
                        switchRows(b_row,t_row); 
 
1498
                        mt.switchRows(b_row,t_row);
 
1499
                }
 
1500
 
 
1501
                //--- now b_row is last non-zero entry from the top
 
1502
                //--- now t_row is first zero entry from the bottom
 
1503
                //--- combine until all entries below diagonal in jTH column = 0
 
1504
 
 
1505
                if ( b_row <= j )
 
1506
                        continue;
 
1507
                        int i;
 
1508
                for( i = j+1; i<= b_row; i++) {
 
1509
                        double f = -mat[i][j] / mat[j][j];
 
1510
                        combineRows(i,f,j);
 
1511
                        mt.combineRows(i,f,j);
 
1512
                }
 
1513
        }
 
1514
        return true;
 
1515
}
 
1516
*/