1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/submat.h,v 1.46 2017/01/12 14:43:54 masarati Exp $ */
2 /*
3 * MBDyn (C) is a multibody analysis code.
4 * http://www.mbdyn.org
5 *
6 * Copyright (C) 1996-2017
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 /* Sottomatrici */
33
34
35 #ifndef SUBMAT_H
36 #define SUBMAT_H
37
38
39 #include <myassert.h>
40 #include <except.h>
41
42
43 /* include del programma */
44 /* si assume che solman.h includa piu' o meno direttamnte f2c.h,
45 * che contiene le dichiarazioni dei tipi derivati dal fortran. */
46
47 #include <solman.h>
48 #include <fullmh.h>
49 #include <matvec3.h>
50 #include <matvec3n.h>
51
52 /* SubMatrixHandler - begin */
53
54 /*
55 Classe virtuale delle sottomatrici.
56 Le SubMatrixHandler sono matrici dotate di vettori di incidenza
57 per righe e colonne.
58 Sono usate per scrivere le sottomatrici di ogni elemento, che poi si
59 sommano con apposite routines alla matrice jacobiana.
60 */
61
62 class SubMatrixHandler : public MatrixHandler {
63 public:
64 /* Errori */
65
66 /*
67 * Errore di ridimensionamento illegale.
68 */
69 class ErrResize : public MBDynErrBase {
70 public:
ErrResize(MBDYN_EXCEPT_ARGS_DECL)71 ErrResize(MBDYN_EXCEPT_ARGS_DECL) : MBDynErrBase(MBDYN_EXCEPT_ARGS_PASSTHRU) {};
72 };
73
74 /* Costruttori */
75
76 /*
77 * Distruttore virtuale.
78 * Necessario per la corretta distruzione degli oggetti derivati
79 */
80 virtual ~SubMatrixHandler(void);
81
82
83 /* Metodi di servizio */
84
85 #ifdef DEBUG
86 /*
87 * Routine di verifica della validita' dell'oggetto.
88 * Usata per il debug.
89 */
90 virtual void IsValid(void) const = 0;
91 #endif /* DEBUG */
92
93
94 /* Inizializzazione */
95
96
97 /*
98 * Ridimensiona la matrice.
99 * Nota: nell'implementazione corrente le dimensioni possono essere solamente
100 * inferiori alle dimensioni massime con cui e' stata dimensionata.
101 */
102 virtual void Resize(integer, integer) = 0;
103
104 /*
105 * Ridimensiona ed inizializza.
106 * Combina le due funzioni precedenti in una chiamata.
107 */
108 virtual void ResizeReset(integer, integer) = 0;
109
110 /* Gestione dei vettori di incidenza */
111
112 /*
113 * Scrive l'indice di riga.
114 */
115 virtual void PutRowIndex(integer, integer) = 0;
116
117 /*
118 * Scrive l'indice di colonna.
119 */
120 virtual void PutColIndex(integer, integer) = 0;
121
122 /*
123 * Ottiene l'indice di riga.
124 */
125 virtual integer iGetRowIndex(integer) const = 0;
126
127 /*
128 * Ottiene l'indice di colonna.
129 */
130 virtual integer iGetColIndex(integer) const = 0;
131
132 /* Funzioni di interazione con le matrici */
133
134 /*
135 * Si somma ad una matrice.
136 * Nota: le dimensioni devono essere compatibili.
137 */
138 virtual MatrixHandler& AddTo(MatrixHandler& MH) const = 0;
139
140 /*
141 * Si somma ad una matrice, trasposta.
142 * Nota: le dimensioni devono essere compatibili.
143 */
144 virtual MatrixHandler& AddToT(MatrixHandler& MH) const = 0;
145
146 /*
147 * Si sottrae da una matrice.
148 * Nota: le dimensioni devono essere compatibili.
149 */
150 virtual MatrixHandler& SubFrom(MatrixHandler& MH) const = 0;
151
152 /*
153 * Si sottrae da una matrice, trasposta.
154 * Nota: le dimensioni devono essere compatibili.
155 */
156 virtual MatrixHandler& SubFromT(MatrixHandler& MH) const = 0;
157 };
158
159 /* SubMatrixHandler - end */
160
161
162 /* FullSubMatrixHandler */
163
164 /*
165 * Sottomatrice piena. E' costituita da un vettore di interi, piRow, che
166 * contiene i due vettori di incidenza, e da un vettore di reali, pdMat,
167 * che contiene la matrice.
168 * Il vettore di incidenza delle righe e' piRow, mentre quello delle colonne
169 * e' piCol = piRow+iNumRows. Le dimensioni possono essere cambiate con
170 * Resize(), con i vincoli che la somma di righe e colonne non deve eccedere
171 * la lunghezza del vettore di interi, ed il loro prodotto non deve eccedere
172 * la lunghezza del vettore di reali.
173 */
174
175 class FullSubMatrixHandler :
176 public SubMatrixHandler, public FullMatrixHandler {
177 friend std::ostream&
178 operator << (std::ostream& out, const FullSubMatrixHandler& m);
179
180 friend class NaiveMatrixHandler;
181 friend class NaivePermMatrixHandler;
182
183
184 protected:
185 /* Dimensione totale del vettore di incidenza */
186 integer iVecSize;
187 /* Puntatore al vettore di incidenza delle righe.
188 * Nota: coincide con il puntatore al vettore di incidenza */
189 integer* piRowm1;
190 /* Puntatore al vettore di incidenza delle colonne */
191 integer* piColm1;
192
193 private:
194 FullSubMatrixHandler(const FullSubMatrixHandler&);
195
196 public:
197
198 /* Costruttori */
199
200 /*
201 * Riceve i vettori per gli indici ed i coefficienti, con relative
202 * dimensioni
203 * @param iIntSize dimensione del vettore degli indici
204 * @param iDoubleSize dimensione del vettore dei coefficienti
205 * @param piTmpVec puntatore al vettore degli indici
206 * @param pdTmpMat puntatore al vettore dei coefficienti
207 */
208 FullSubMatrixHandler(integer iIntSize, integer* piTmpVec,
209 integer iDoubleSize, doublereal* pdTmpMat,
210 integer iMaxCols, doublereal **ppdCols);
211
212 FullSubMatrixHandler(integer iNR, integer iNC = 0);
213
214 /*
215 * Distruttore banale.
216 * Nota: l'elemento non possiede memoria e quindi non ne dealloca.
217 */
218 virtual ~FullSubMatrixHandler(void);
219
220 /* Metodi di servizio */
221
222 #ifdef DEBUG
223 /*
224 * Routine di verifica della validita' dell'oggetto.
225 * Usata per il debug.
226 */
227 virtual void IsValid(void) const;
228 #endif /* DEBUG */
229
230 /*
231 * Numero di righe della sottomatrice
232 */
iGetNumRows(void)233 integer iGetNumRows(void) const {
234 return iNumRows;
235 };
236
237 /*
238 * Numero di colonne della sottomatrice
239 */
iGetNumCols(void)240 integer iGetNumCols(void) const {
241 return iNumCols;
242 };
243
244 /* Metodi di inizializzazione */
245
246 /*
247 * Inizializza la porzione utilizzata con il valore desiderato
248 */
249 void Reset(void);
250
251 /*
252 * Modifica le dimensioni correnti
253 */
254 void Resize(integer iNewRow, integer iNewCol);
255
256 /* Ridimensiona ed inizializza. */
257 virtual void ResizeReset(integer, integer);
258
259 /*
260 * Collega la matrice Full alla memoria che gli viene passata
261 * in ingresso
262 */
263 void Attach(int iRows, int iCols, integer* piTmpIndx);
264
265 /* Gestione dei coefficienti */
266
267 /* Inserisce un coefficiente */
268
269 /*
270 * NOTE:
271 *
272 * these functions need be redefined here because their
273 * inheritance is ambiguous, and we want to use the
274 * (very efficient) FullMatrixHandler version instead
275 * of the (less efficient) MatrixHandler version, which
276 * is a wrapper for the () operator.
277 *
278 * However, if compiled with appropriate optimization,
279 * they are absolutely equivalent to directly accessing
280 * the matrix, as the #if 0'ed code does.
281 */
282 inline void
283 PutCoef(integer iRow, integer iCol, const doublereal& dCoef);
284
285 /* Incrementa un coefficiente - se non esiste lo crea */
286 inline void
287 IncCoef(integer iRow, integer iCol, const doublereal& dCoef);
288
289 /* Decrementa un coefficiente - se non esiste lo crea */
290 inline void
291 DecCoef(integer iRow, integer iCol, const doublereal& dCoef);
292
293 /* Restituisce un coefficiente - zero se non e' definito */
294 inline const doublereal&
295 dGetCoef(integer iRow, integer iCol) const;
296
297 // FIXME: disambiguate operator()
298 inline const doublereal&
299 operator () (integer iRow, integer iCol) const;
300
301 inline doublereal&
302 operator () (integer iRow, integer iCol);
303 // end of FIXME: disambiguate operator()
304
305 /* Gestione degli indici */
306
307 /*
308 * Scrive un indice di riga
309 */
310 inline void
PutRowIndex(integer iSubRow,integer iRow)311 PutRowIndex(integer iSubRow, integer iRow) {
312 #ifdef DEBUG
313 IsValid();
314 #endif /* DEBUG */
315
316 ASSERT((iSubRow > 0) && (iSubRow <= iNumRows));
317
318 piRowm1[iSubRow] = iRow;
319 };
320
321 /*
322 * Scrive un indice di colonna
323 */
324 inline void
PutColIndex(integer iSubCol,integer iCol)325 PutColIndex(integer iSubCol, integer iCol) {
326 #ifdef DEBUG
327 IsValid();
328 #endif /* DEBUG */
329
330 ASSERT((iSubCol > 0) && (iSubCol <= iNumCols));
331
332 piColm1[iSubCol] = iCol;
333 };
334
335 /*
336 * Legge un indice di riga
337 */
338 inline integer
iGetRowIndex(integer iSubRow)339 iGetRowIndex(integer iSubRow) const {
340 #ifdef DEBUG
341 IsValid();
342 #endif /* DEBUG */
343
344 ASSERT((iSubRow > 0) && (iSubRow <= iNumRows));
345
346 return piRowm1[iSubRow];
347 };
348
349 /*
350 * Legge un indice di colonna
351 */
352 inline integer
iGetColIndex(integer iSubCol)353 iGetColIndex(integer iSubCol) const {
354 #ifdef DEBUG
355 IsValid();
356 #endif /* DEBUG */
357
358 ASSERT((iSubCol > 0) && (iSubCol <= iNumCols));
359
360 return piColm1[iSubCol];
361 };
362
363 /*
364 * Somma un vettore di tipo Vec3 in una data posizione.
365 * Nota: si assume che nella sottomatrice vi sia spazio
366 * per il vettore 3.
367 * Nota: gli indici sono a base 1, in stile FORTRAN.
368 * @param iRow indice di riga della sottomatrice da cui iniziare
369 * @param iCol indice di colonna della sottomatrice da cui iniziare
370 * @param v Vec3 da sommare
371 */
372 void
373 Add(integer iRow, integer iCol, const Vec3& v);
374
375 /*
376 * Sottrae un vettore di tipo Vec3 in una data posizione.
377 * Nota: si assume che nella sottomatrice vi sia spazio
378 * per il vettore 3.
379 * Nota: gli indici sono a base 1, in stile FORTRAN.
380 * @param iRow indice di riga della sottomatrice da cui iniziare
381 * @param iCol indice di colonna della sottomatrice da cui iniziare
382 * @param v Vec3 da sommare
383 */
384 void
385 Sub(integer iRow, integer iCol, const Vec3& v);
386
387 /*
388 * Scrive un vettore di tipo Vec3 in una data posizione.
389 * Nota: si assume che nella sottomatrice vi sia spazio
390 * per il vettore 3.
391 * Nota: gli indici sono a base 1, in stile FORTRAN.
392 * @param iRow indice di riga della sottomatrice da cui iniziare
393 * @param iCol indice di colonna della sottomatrice da cui iniziare
394 * @param v Vec3 da sommare
395 */
396 void
397 Put(integer iRow, integer iCol, const Vec3& v);
398
399 /*
400 * Somma un vettore di tipo Vec3 trasposto in una data posizione.
401 * Nota: si assume che nella sottomatrice vi sia spazio
402 * per il vettore 3 trasposto (riga).
403 * Nota: gli indici sono a base 1, in stile FORTRAN.
404 * @param iRow indice di riga della sottomatrice da cui iniziare
405 * @param iCol indice di colonna della sottomatrice da cui iniziare
406 * @param v Vec3 da sommare
407 */
408 void
409 AddT(integer iRow, integer iCol, const Vec3& v);
410
411 /*
412 * Sottrae un vettore di tipo Vec3 trasposto in una data posizione.
413 * Nota: si assume che nella sottomatrice vi sia spazio
414 * per il vettore 3 trasposto (riga).
415 * Nota: gli indici sono a base 1, in stile FORTRAN.
416 * @param iRow indice di riga della sottomatrice da cui iniziare
417 * @param iCol indice di colonna della sottomatrice da cui iniziare
418 * @param v Vec3 da sommare
419 */
420 void
421 SubT(integer iRow, integer iCol, const Vec3& v);
422
423 /*
424 * Scrive un vettore di tipo Vec3 trasposto in una data posizione.
425 * Nota: si assume che nella sottomatrice vi sia spazio
426 * per il vettore 3 trasposto (riga).
427 * Nota: gli indici sono a base 1, in stile FORTRAN.
428 * @param iRow indice di riga della sottomatrice da cui iniziare
429 * @param iCol indice di colonna della sottomatrice da cui iniziare
430 * @param v Vec3 da sommare
431 */
432 void
433 PutT(integer iRow, integer iCol, const Vec3& v);
434
435 #if 0 /* FIXME: replace original? */
436 /*
437 * Somma un vettore di tipo Vec3 in una data posizione in diagonale.
438 * Nota: si assume che nella sottomatrice vi sia spazio
439 * per il vettore 3.
440 * Nota: gli indici sono a base 1, in stile FORTRAN.
441 * @param iRow indice di riga della sottomatrice da cui iniziare
442 * @param iCol indice di colonna della sottomatrice da cui iniziare
443 * @param v Vec3 da sommare
444 */
445 void
446 AddDiag(integer iRow, integer iCol, const Vec3& v);
447
448 /*
449 * Sottrae un vettore di tipo Vec3 in una data posizione in diagonale.
450 * Nota: si assume che nella sottomatrice vi sia spazio
451 * per il vettore 3.
452 * Nota: gli indici sono a base 1, in stile FORTRAN.
453 * @param iRow indice di riga della sottomatrice da cui iniziare
454 * @param iCol indice di colonna della sottomatrice da cui iniziare
455 * @param v Vec3 da sommare
456 */
457 void
458 SubDiag(integer iRow, integer iCol, const Vec3& v);
459
460 /*
461 * Scrive un vettore di tipo Vec3 in una data posizione in diagonale.
462 * Nota: si assume che nella sottomatrice vi sia spazio
463 * per il vettore 3.
464 * Nota: gli indici sono a base 1, in stile FORTRAN.
465 * @param iRow indice di riga della sottomatrice da cui iniziare
466 * @param iCol indice di colonna della sottomatrice da cui iniziare
467 * @param v Vec3 da sommare
468 */
469 void
470 PutDiag(integer iRow, integer iCol, const Vec3& v);
471
472 /*
473 * Somma un vettore di tipo Vec3 in una data posizione [ v x ].
474 * Nota: si assume che nella sottomatrice vi sia spazio
475 * per il vettore 3.
476 * Nota: gli indici sono a base 1, in stile FORTRAN.
477 * @param iRow indice di riga della sottomatrice da cui iniziare
478 * @param iCol indice di colonna della sottomatrice da cui iniziare
479 * @param v Vec3 da sommare
480 */
481 void
482 AddCross(integer iRow, integer iCol, const Vec3& v);
483
484 /*
485 * Sottrae un vettore di tipo Vec3 in una data posizione [ v x ].
486 * Nota: si assume che nella sottomatrice vi sia spazio
487 * per il vettore 3.
488 * Nota: gli indici sono a base 1, in stile FORTRAN.
489 * @param iRow indice di riga della sottomatrice da cui iniziare
490 * @param iCol indice di colonna della sottomatrice da cui iniziare
491 * @param v Vec3 da sommare
492 */
493 void
494 SubCross(integer iRow, integer iCol, const Vec3& v);
495
496 /*
497 * Scrive un vettore di tipo Vec3 in una data posizione [ v x ].
498 * Nota: si assume che nella sottomatrice vi sia spazio
499 * per il vettore 3.
500 * Nota: gli indici sono a base 1, in stile FORTRAN.
501 * @param iRow indice di riga della sottomatrice da cui iniziare
502 * @param iCol indice di colonna della sottomatrice da cui iniziare
503 * @param v Vec3 da sommare
504 */
505 void
506 PutCross(integer iRow, integer iCol, const Vec3& v);
507 #endif
508
509 /*
510 * Somma una matrice di tipo Mat3x3 in una data posizione.
511 * Nota: si assume che nella sottomatrice vi sia spazio
512 * per la matrice 3x3.
513 * Nota: gli indici sono a base 1, in stile FORTRAN.
514 * @param iRow indice di riga della sottomatrice da cui iniziare
515 * @param iCol indice di colonna della sottomatrice da cui iniziare
516 * @param m Mat3x3 da sommare
517 */
518 void
519 Add(integer iRow, integer iCol, const Mat3x3& m);
520 void
521 AddT(integer iRow, integer iCol, const Mat3x3& m);
522
523 /*
524 * Sottrae una matrice di tipo Mat3x3 da una data posizione.
525 * Nota: si assume che nella sottomatrice vi sia spazio
526 * per la matrice 3x3.
527 * Nota: gli indici sono a base 1, in stile FORTRAN.
528 * @param iRow indice di riga della sottomatrice da cui iniziare
529 * @param iCol indice di colonna della sottomatrice da cui iniziare
530 * @param m Mat3x3 da sottrarre
531 */
532 void
533 Sub(integer iRow, integer iCol, const Mat3x3& m);
534 void
535 SubT(integer iRow, integer iCol, const Mat3x3& m);
536
537 /*
538 * Scrive una matrice di tipo Mat3x3 in una data posizione.
539 * Nota: si assume che nella sottomatrice vi sia spazio
540 * per la matrice 3x3.
541 * Nota: gli indici sono a base 1, in stile FORTRAN.
542 * @param iRow indice di riga della sottomatrice da cui iniziare
543 * @param iCol indice di colonna della sottomatrice da cui iniziare
544 * @param m Mat3x3 da scrivere
545 */
546 void
547 Put(integer iRow, integer iCol, const Mat3x3& m);
548 void
549 PutT(integer iRow, integer iCol, const Mat3x3& m);
550
551 /*
552 * Somma una matrice di tipo Mat3xN in una data posizione.
553 * Nota: si assume che nella sottomatrice vi sia spazio
554 * per la matrice 3x3.
555 * Nota: gli indici sono a base 1, in stile FORTRAN.
556 * @param iRow indice di riga della sottomatrice da cui iniziare
557 * @param iCol indice di colonna della sottomatrice da cui iniziare
558 * @param m Mat3xN
559 */
560 void
561 Add(integer iRow, integer iCol, const Mat3xN& m);
562 void
563 AddT(integer iRow, integer iCol, const Mat3xN& m);
564
565 /*
566 * Sottrae una matrice di tipo Mat3xN da una data posizione.
567 * Nota: si assume che nella sottomatrice vi sia spazio
568 * per la matrice 3x3.
569 * Nota: gli indici sono a base 1, in stile FORTRAN.
570 * @param iRow indice di riga della sottomatrice da cui iniziare
571 * @param iCol indice di colonna della sottomatrice da cui iniziare
572 * @param m Mat3xN
573 */
574 void
575 Sub(integer iRow, integer iCol, const Mat3xN& m);
576 void
577 SubT(integer iRow, integer iCol, const Mat3xN& m);
578
579 /*
580 * Scrive una matrice di tipo Mat3xN in una data posizione.
581 * Nota: si assume che nella sottomatrice vi sia spazio
582 * per la matrice 3x3.
583 * Nota: gli indici sono a base 1, in stile FORTRAN.
584 * @param iRow indice di riga della sottomatrice da cui iniziare
585 * @param iCol indice di colonna della sottomatrice da cui iniziare
586 * @param m Mat3xN
587 */
588 void
589 Put(integer iRow, integer iCol, const Mat3xN& m);
590 void
591 PutT(integer iRow, integer iCol, const Mat3xN& m);
592
593 /* come sopra, ma per matrici Nx3 **/
594 void Add(integer iRow, integer iCol, const MatNx3& m);
595 void Sub(integer iRow, integer iCol, const MatNx3& m);
596 void Put(integer iRow, integer iCol, const MatNx3& m);
597
598 void PutDiag(integer iFirstRow, integer iFirstCol, const Vec3& v);
599 void PutDiag(integer iFirstRow, integer iFirstCol, const doublereal& v);
600 void PutCross(integer iFirstRow, integer iFirstCol, const Vec3& v);
601
602 void Add(integer iRow, integer iCol,
603 const FullMatrixHandler & source);
604 void Sub(integer iRow, integer iCol,
605 const FullMatrixHandler & source);
606 void Put(integer iRow, integer iCol,
607 const FullMatrixHandler & source);
608 void Add(integer iRow, integer iCol,
609 const FullMatrixHandler & source, const doublereal dCoef);
610 void Sub(integer iRow, integer iCol,
611 const FullMatrixHandler & source, const doublereal dCoef);
612 void Put(integer iRow, integer iCol,
613 const FullMatrixHandler & source, const doublereal dCoef);
614 void AddT(integer iRow, integer iCol,
615 const FullMatrixHandler & source);
616 void SubT(integer iRow, integer iCol,
617 const FullMatrixHandler & source);
618 void PutT(integer iRow, integer iCol,
619 const FullMatrixHandler & source);
620 void AddT(integer iRow, integer iCol,
621 const FullMatrixHandler & source, const doublereal dCoef);
622 void SubT(integer iRow, integer iCol,
623 const FullMatrixHandler & source, const doublereal dCoef);
624 void PutT(integer iRow, integer iCol,
625 const FullMatrixHandler & source, const doublereal dCoef);
626
627 /* Interazione con le matrici */
628
629 /*
630 * Somma la matrice ad un matrix handler usando i metodi generici
631 */
632 MatrixHandler& AddTo(MatrixHandler& MH) const;
633
634 /*
635 * Somma la matrice, trasposta, ad un matrix handler usando i metodi generici
636 */
637 MatrixHandler& AddToT(MatrixHandler& MH) const;
638
639 /*
640 * Somma la matrice ad un FullMatrixHandler
641 */
642 MatrixHandler& AddTo(FullMatrixHandler& MH) const;
643
644 /*
645 * Somma la matrice, trasposta, ad un FullMatrixHandler
646 */
647 MatrixHandler& AddToT(FullMatrixHandler& MH) const;
648
649 /*
650 * Sottrae la matrice da un matrix handler usando i metodi generici
651 */
652 MatrixHandler& SubFrom(MatrixHandler& MH) const;
653
654 /*
655 * Sottrae la matrice, trasposta, da un matrix handler usando i metodi generici
656 */
657 MatrixHandler& SubFromT(MatrixHandler& MH) const;
658
659 /*
660 * Sottrae la matrice da un FullMatrixHandler
661 */
662 MatrixHandler& SubFrom(FullMatrixHandler& MH) const;
663
664 /*
665 * Sottrae la matrice, trasposta, da un FullMatrixHandler
666 */
667 MatrixHandler& SubFromT(FullMatrixHandler& MH) const;
668 };
669
670 /* Inserisce un coefficiente */
671 inline void
PutCoef(integer iRow,integer iCol,const doublereal & dCoef)672 FullSubMatrixHandler::PutCoef(integer iRow, integer iCol,
673 const doublereal& dCoef)
674 {
675 #if 0
676 ppdColsm1[iCol][iRow] = dCoef;
677 #endif
678 FullMatrixHandler::PutCoef(iRow, iCol, dCoef);
679 }
680
681 /* Incrementa un coefficiente - se non esiste lo crea */
682 inline void
IncCoef(integer iRow,integer iCol,const doublereal & dCoef)683 FullSubMatrixHandler::IncCoef(integer iRow, integer iCol,
684 const doublereal& dCoef)
685 {
686 #if 0
687 ppdColsm1[iCol][iRow] += dCoef;
688 #endif
689 FullMatrixHandler::IncCoef(iRow, iCol, dCoef);
690 }
691
692 /* Decrementa un coefficiente - se non esiste lo crea */
693 inline void
DecCoef(integer iRow,integer iCol,const doublereal & dCoef)694 FullSubMatrixHandler::DecCoef(integer iRow, integer iCol,
695 const doublereal& dCoef)
696 {
697 #if 0
698 ppdColsm1[iCol][iRow] -= dCoef;
699 #endif
700 FullMatrixHandler::DecCoef(iRow, iCol, dCoef);
701 }
702
703 /* Restituisce un coefficiente - zero se non e' definito */
704 inline const doublereal&
dGetCoef(integer iRow,integer iCol)705 FullSubMatrixHandler::dGetCoef(integer iRow, integer iCol) const
706 {
707 #if 0
708 return ppdColsm1[iCol][iRow];
709 #endif
710 return FullMatrixHandler::dGetCoef(iRow, iCol);
711 }
712
713 // FIXME: disambiguate operator()
714 inline const doublereal&
operator()715 FullSubMatrixHandler::operator () (integer iRow, integer iCol) const
716 {
717 return FullMatrixHandler::operator()(iRow, iCol);
718 }
719
720 inline doublereal&
operator()721 FullSubMatrixHandler::operator () (integer iRow, integer iCol)
722 {
723 return FullMatrixHandler::operator()(iRow, iCol);
724 }
725 // end of FIXME: disambiguate operator()
726
727 /* FullSubMatrixHandler - end */
728
729
730 /* SparseSubMatrixHandler */
731
732 /*
733 * Gestore di sottomatrici sparse, piuttosto rozzo, va usato con cautela.
734 * E' formato da due vettori di interi e da uno di reali, tutti della stessa
735 * lunghezza.
736 * Ad ogni indice del sotto-vettore corrisponde un coefficiente con i suoi
737 * due indici nella matrice completa.
738 * La scrittura non e' strutturata, per cui l'utilizzatore deve badare a non
739 * lasciare vuoti e a non ripetere i coefficienti (in realta' non succede
740 * nulla se la si usa per assemblare, solo uno stesso coefficiente puo' essere
741 * dato da piu' contributi).
742 */
743
744 class SparseSubMatrixHandler : public SubMatrixHandler {
745 friend class SparseMatrixHandler;
746 friend class FullMatrixHandler;
747 friend class NaiveMatrixHandler;
748 friend class NaivePermMatrixHandler;
749
750 public:
751 /* Errori */
752
753 class ErrResize : MBDynErrBase {
754 public:
ErrResize(MBDYN_EXCEPT_ARGS_DECL)755 ErrResize(MBDYN_EXCEPT_ARGS_DECL) : MBDynErrBase(MBDYN_EXCEPT_ARGS_PASSTHRU) {};
756 };
757
758 private:
759 bool bOwnsMemory;
760 /* Dimensioni dell'array degli indici */
761 integer iIntSize;
762 /* Dimensioni dell'array dei coefficienti */
763 integer iDoubleSize;
764 /* Numero di entries definite */
765 integer iNumItems;
766 /* Puntatore all'array degli indici di riga.
767 * Coincide con il puntatore all'array degli inidici, che e' unico */
768 integer* piRowm1;
769 /* Puntatore all'array degli indici di colonna */
770 integer* piColm1;
771 /* Puntatore all'array dei coefficienti */
772 doublereal* pdMatm1;
773
774 private:
775 SparseSubMatrixHandler(const SparseSubMatrixHandler&);
776
777 public:
778 /* Costruttori */
779
780
781 /* Costruttore.
782 * @param iTmpInt dimensione dell'array degli indici
783 * @param iTmpDouble dimensione dell'array dei coefficienti
784 * @param piTmpIndex puntatore all'array degli indici
785 * @param pdTmpMat puntatore all'array dei coefficienti
786 */
787 SparseSubMatrixHandler(integer iTmpInt, integer* piTmpIndex,
788 integer iTmpDouble, doublereal* pdTmpMat);
789
790 SparseSubMatrixHandler(integer iTmpInt);
791
792 /* Distruttore banale.
793 * Nota: dato che la classe non possiede la memoria,
794 * non ne deve deallocare
795 */
796 virtual ~SparseSubMatrixHandler(void);
797
798 /* Metodi di servizio */
799
800 #ifdef DEBUG
801 /*
802 * Routine di verifica della validita' dell'oggetto.
803 * Usata per il debug.
804 */
805 virtual void IsValid(void) const;
806 #endif /* DEBUG */
807
808 /*
809 * Numero di righe della sottomatrice.
810 * Nota: rappresenta il numero totale di entries della sottomatrice.
811 */
iGetNumRows(void)812 integer iGetNumRows(void) const {
813 return iNumItems;
814 };
815
816 /*
817 * Numero di colonne della sottomatrice.
818 * Nota: e' sempre 1, ovvero la matrice e' interpretata
819 * come un vettore.
820 */
iGetNumCols(void)821 integer iGetNumCols(void) const {
822 return 1;
823 };
824
825 /* Metodi di inizializzazione */
826
827 /*
828 * Ridimensiona la matrice.
829 * Nota: solo il primo argomento viene considerato,
830 * e rappresenta il numero totale di entries.
831 * Questo metodo deve essere chiamato prima di qualsiasi
832 * operazione sulla matrice.
833 */
834 void Resize(integer iNewRow, integer iNewCol);
835
836 /*
837 * Ridimensiona ed inizializza.
838 * Unione dei due metodi precedenti
839 */
840 void ResizeReset(integer iNewRow, integer iNewCol);
841
842 /* Azzera */
843 void Reset(void);
844
845 /*
846 * Collega la matrice sparsa alla memoria che gli viene passata
847 * in ingresso
848 */
849 void Attach(int iNumEntr, doublereal* pdTmpMat, integer* piTmpIndx);
850
851 /* Gestione dei coefficienti */
852
853 /*
854 * Scrive un coefficiente in base ai sottoindici.
855 */
856 inline void
PutCoef(integer iSubIt,integer iDmy,const doublereal & dCoef)857 PutCoef(integer iSubIt, integer iDmy, const doublereal& dCoef) {
858 #ifdef DEBUG
859 IsValid();
860 #endif /* DEBUG */
861
862 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
863
864 pdMatm1[iSubIt] = dCoef;
865 };
866
867 /*
868 * Incrementa un coefficiente in base ai sottoindici.
869 */
870 inline void
IncCoef(integer iSubIt,integer iDmy,const doublereal & dCoef)871 IncCoef(integer iSubIt, integer iDmy, const doublereal& dCoef) {
872 #ifdef DEBUG
873 IsValid();
874 #endif /* DEBUG */
875
876 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
877 pdMatm1[iSubIt] += dCoef;
878 };
879
880 /*
881 * Decrementa un coefficiente in base ai sottoindici.
882 */
883 inline void
DecCoef(integer iSubIt,integer iDmy,const doublereal & dCoef)884 DecCoef(integer iSubIt, integer iDmy, const doublereal& dCoef) {
885 #ifdef DEBUG
886 IsValid();
887 #endif /* DEBUG */
888
889 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
890 pdMatm1[iSubIt] -= dCoef;
891 };
892
893 /*
894 * Ottiene un coefficiente in base ai sottoindici.
895 */
896 inline const doublereal&
dGetCoef(integer iSubIt,integer iDmy)897 dGetCoef(integer iSubIt, integer iDmy) const {
898 #ifdef DEBUG
899 IsValid();
900 #endif /* DEBUG */
901
902 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
903
904 return pdMatm1[iSubIt];
905 };
906
907 /*
908 * Ottiene un coefficiente in base ai sottoindici.
909 */
910 inline const doublereal&
operator()911 operator () (integer iSubIt, integer iDmy) const {
912 #ifdef DEBUG
913 IsValid();
914 #endif /* DEBUG */
915
916 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
917
918 return pdMatm1[iSubIt];
919 };
920
921 /*
922 * Ottiene un coefficiente in base ai sottoindici.
923 */
924 inline doublereal&
operator()925 operator () (integer iSubIt, integer iDmy) {
926 #ifdef DEBUG
927 IsValid();
928 #endif /* DEBUG */
929
930 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
931
932 return pdMatm1[iSubIt];
933 };
934
935 /*
936 * Scrive un indice di riga
937 */
938 inline void
PutRowIndex(integer iSubIt,integer iRow)939 PutRowIndex(integer iSubIt, integer iRow) {
940 #ifdef DEBUG
941 IsValid();
942 #endif /* DEBUG */
943
944 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
945 piRowm1[iSubIt] = iRow;
946 };
947
948 /*
949 * Scrive un indice di colonna
950 */
951 inline void
PutColIndex(integer iSubIt,integer iCol)952 PutColIndex(integer iSubIt, integer iCol) {
953 #ifdef DEBUG
954 IsValid();
955 #endif /* DEBUG */
956
957 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
958 piColm1[iSubIt] = iCol;
959 };
960
961 /*
962 * Ottiene un indice di riga
963 */
964 inline integer
iGetRowIndex(integer iSubIt)965 iGetRowIndex(integer iSubIt) const {
966 #ifdef DEBUG
967 IsValid();
968 #endif /* DEBUG */
969
970 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
971
972 return piRowm1[iSubIt];
973 };
974
975 /*
976 * Ottiene un indice di colonna
977 */
978 inline integer
iGetColIndex(integer iSubIt)979 iGetColIndex(integer iSubIt) const {
980 #ifdef DEBUG
981 IsValid();
982 #endif /* DEBUG */
983
984 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
985
986 return piColm1[iSubIt];
987 };
988
989 /*
990 * Scrive un'entry completa.
991 * @param iSubIt sottoindice (numero della entry)
992 * @param iRow indice di riga
993 * @param iCol indice di colonna
994 * @param dCoef coefficiente
995 */
996 inline void
PutItem(integer iSubIt,integer iRow,integer iCol,const doublereal & dCoef)997 PutItem(integer iSubIt, integer iRow, integer iCol,
998 const doublereal& dCoef) {
999 #ifdef DEBUG
1000 IsValid();
1001 #endif /* DEBUG */
1002
1003 ASSERT((iSubIt > 0) && (iSubIt <= iNumItems));
1004 ASSERT(iRow > 0);
1005 ASSERT(iCol > 0);
1006
1007 pdMatm1[iSubIt] = dCoef;
1008 piRowm1[iSubIt] = iRow;
1009 piColm1[iSubIt] = iCol;
1010 };
1011
1012 /*
1013 * Scrive una matrice diagonale nella posizione assegnata.
1014 * @param iSubIt sottoindice iniziale (numero della prima entry)
1015 * @param iFirstRow indice della prima riga della matrice completa
1016 * @param iFirstCol indice della prima colonna della matrice completa
1017 * @param v vettore diagonale della matrice
1018 */
1019 void
1020 PutDiag(integer iSubIt, integer iFirstRow, integer iFirstCol,
1021 const Vec3& v);
1022
1023 /*
1024 * Scrive una matrice diagonale nella posizione assegnata.
1025 * @param iSubIt sottoindice iniziale (numero della prima entry)
1026 * @param iFirstRow indice della prima riga della matrice completa
1027 * @param iFirstCol indice della prima colonna della matrice completa
1028 * @param d coefficiente della diagonale della matrice
1029 */
1030 void
1031 PutDiag(integer iSubIt, integer iFirstRow, integer iFirstCol,
1032 const doublereal& d);
1033
1034 /*
1035 * Scrive una matrice prodotto vettore nella posizione assegnata.
1036 * @param iSubIt sottoindice iniziale (numero della prima entry)
1037 * @param iFirstRow indice della prima riga della matrice completa
1038 * @param iFirstCol indice della prima colonna della matrice completa
1039 * @param v vettore da cui viene calcolata la matrice prodotto
1040 * vettore
1041 */
1042 void
1043 PutCross(integer iSubIt, integer iFirstRow, integer iFirstCol,
1044 const Vec3& v);
1045
1046 /*
1047 * Scrive una Mat3x3 nella posizione assegnata.
1048 * @param iSubIt sottoindice iniziale (numero della prima entry)
1049 * @param iFirstRow indice della prima riga della matrice completa
1050 * @param iFirstCol indice della prima colonna della matrice completa
1051 * @param m matrice da inserire
1052 */
1053 void
1054 PutMat3x3(integer iSubIt, integer iFirstRow, integer iFirstCol,
1055 const Mat3x3& m);
1056
1057 /* Interazione con le matrici */
1058
1059 /*
1060 * Somma la matrice ad un matrix handler usando i metodi generici
1061 */
1062 MatrixHandler& AddTo(MatrixHandler& MH) const;
1063
1064 /*
1065 * Somma la matrice, trasposta, ad un matrix handler usando i metodi generici
1066 */
1067 MatrixHandler& AddToT(MatrixHandler& MH) const;
1068
1069 /*
1070 * Somma la matrice ad un FullMatrixHandler
1071 */
1072 MatrixHandler& AddTo(FullMatrixHandler& MH) const;
1073
1074 /*
1075 * Somma la matrice, trasposta, ad un FullMatrixHandler
1076 */
1077 MatrixHandler& AddToT(FullMatrixHandler& MH) const;
1078
1079 /*
1080 * Sottrae la matrice da un matrix handler usando i metodi generici
1081 */
1082 MatrixHandler& SubFrom(MatrixHandler& MH) const;
1083
1084 /*
1085 * Sottrae la matrice, trasposta, da un matrix handler usando i metodi generici
1086 */
1087 MatrixHandler& SubFromT(MatrixHandler& MH) const;
1088
1089 /*
1090 * Sottrae la matrice da un FullMatrixHandler
1091 */
1092 MatrixHandler& SubFrom(FullMatrixHandler& MH) const;
1093
1094 /*
1095 * Sottrae la matrice, trasposta, da un FullMatrixHandler
1096 */
1097 MatrixHandler& SubFromT(FullMatrixHandler& MH) const;
1098 };
1099
1100 /* SparseSubMatrixHandler - end */
1101
1102
1103 /* VariableSubMatrixHandler - begin */
1104
1105 /*
1106 * Matrice che puo' diventare via via una sottomatrice piena o una sottomatrice
1107 * sparsa, condividendo la memoria.
1108 * Viene passata agli elementi che, a seconda della loro convenienza,
1109 * la configurano nel modo piu' opportuno.
1110 * Quindi, con metodi opportuni, viene sommata alla matrice completa.
1111 */
1112
1113 class VariableSubMatrixHandler
1114 : public FullSubMatrixHandler, public SparseSubMatrixHandler {
1115 friend class NaiveMatrixHandler;
1116 friend class NaivePermMatrixHandler;
1117 private:
1118 /*
1119 * Stato della matrice.
1120 */
1121 enum { NULLMATRIX, FULL, SPARSE } eStatus;
1122
1123 VariableSubMatrixHandler(const VariableSubMatrixHandler&);
1124
1125 public:
1126 /* Costruttori */
1127
1128 /*
1129 * Costruttore: riceve gli spazi di lavoro con le loro dimensioni
1130 * ed inizializza le matrici piena e sparsa.
1131 * @param iIntSize dimensioni dell'array degli indici
1132 * @param iDoubleSize dimensioni dell'array dei coefficienti
1133 * @param piInt array degli indici
1134 * @param pdDouble array dei coefficienti
1135 */
VariableSubMatrixHandler(integer iIntSize,integer * piInt,integer iDoubleSize,doublereal * pdDouble,integer iMaxRows,integer iMaxCols)1136 VariableSubMatrixHandler(integer iIntSize, integer* piInt,
1137 integer iDoubleSize, doublereal* pdDouble,
1138 integer iMaxRows, integer iMaxCols)
1139 : FullSubMatrixHandler(iMaxRows, iMaxCols),
1140 SparseSubMatrixHandler(iIntSize, piInt, iDoubleSize, pdDouble),
1141 eStatus(NULLMATRIX) {
1142 NO_OP;
1143 };
1144
1145 VariableSubMatrixHandler(integer iMaxRows, integer iMaxCols, integer iNumItems = -1)
FullSubMatrixHandler(iMaxRows,iMaxCols)1146 : FullSubMatrixHandler(iMaxRows, iMaxCols),
1147 SparseSubMatrixHandler(iNumItems >= 0 ? iNumItems : iMaxRows * iMaxCols),
1148 eStatus(NULLMATRIX)
1149 {
1150 NO_OP;
1151 };
1152
1153 /* Metodi di servizio */
1154
1155 /*
1156 * Setta la matrice come vuota.
1157 * Di conseguenza non viene assemblata.
1158 */
SetNullMatrix(void)1159 void SetNullMatrix(void) {
1160 eStatus = NULLMATRIX;
1161 };
1162
1163 /*
1164 * Setta la matrice come piena.
1165 * Ritorna un riferimento a matrice piena, che puo' essere usato
1166 * per le normali operazioni di scrittura delle matrici piene.
1167 */
SetFull(void)1168 FullSubMatrixHandler& SetFull(void) {
1169 eStatus = FULL;
1170 return *dynamic_cast<FullSubMatrixHandler *>(this);
1171 };
1172
1173 /*
1174 * Setta la matrice come sparsa.
1175 * Ritorna un riferimento a matrice sparsa, che puo' essere usato
1176 * per le normali operazioni di scrittura delle matrici sparse.
1177 */
SetSparse(void)1178 SparseSubMatrixHandler& SetSparse(void) {
1179 eStatus = SPARSE;
1180 return *dynamic_cast<SparseSubMatrixHandler *>(this);
1181 };
1182
1183 /*
1184 * Verifica se la matrice e' vuota.
1185 */
bIsNullMatrix(void)1186 bool bIsNullMatrix(void) const {
1187 return (eStatus == NULLMATRIX);
1188 };
1189
1190 /*
1191 * Verifica se la matrice e' piena.
1192 */
bIsFull(void)1193 bool bIsFull(void) const {
1194 return (eStatus == FULL);
1195 };
1196
1197 /*
1198 * Verifica se la matrice e' sparsa.
1199 */
bIsSparse(void)1200 bool bIsSparse(void) const {
1201 return (eStatus == SPARSE);
1202 };
1203
1204 #if 0
1205 /*
1206 * Numero di righe della sottomatrice
1207 */
1208 integer iGetNumRows(void) const {
1209 switch (eStatus) {
1210 case FULL:
1211 return FullSubMatrixHandler::iGetNumRows();
1212
1213 case SPARSE:
1214 return SparseSubMatrixHandler::iGetNumRows();
1215
1216 default:
1217 return 0;
1218 }
1219 };
1220
1221 /*
1222 * Numero di colonne della sottomatrice
1223 */
1224 integer iGetNumCols(void) const {
1225 switch (eStatus) {
1226 case FULL:
1227 return FullSubMatrixHandler::iGetNumCols();
1228
1229 case SPARSE:
1230 return SparseSubMatrixHandler::iGetNumCols();
1231
1232 default:
1233 return 0;
1234 }
1235 };
1236
1237 /*
1238 * Links sparse matrix with already assigned memory
1239 */
1240 void Attach(int iNumEntr, doublereal* pdTmpMat, integer* piTmpIndx) {
1241 SetSparse();
1242 SparseSubMatrixHandler::Attach(iNumEntr, pdTmpMat, piTmpIndx);
1243 };
1244
1245 /*
1246 * Links full matrix with already assigned memory
1247 */
1248 void Attach(int iNumRows, int iNumCols,
1249 doublereal* pdTmpMat, integer* piTmpIndx) {
1250 SetFull();
1251 FullSubMatrixHandler::Attach(iNumRows, iNumCols,
1252 pdTmpMat, piTmpIndx);
1253 };
1254 #endif
1255
1256 /* Interazione con le matrici */
1257
1258 /*
1259 * Si somma ad una matrice completa con metodi generici.
1260 */
AddTo(MatrixHandler & MH)1261 MatrixHandler& AddTo(MatrixHandler& MH) const {
1262 switch (eStatus) {
1263 case FULL:
1264 return FullSubMatrixHandler::AddTo(MH);
1265
1266 case SPARSE:
1267 return SparseSubMatrixHandler::AddTo(MH);
1268
1269 default:
1270 return MH;
1271 }
1272 };
1273
1274 /*
1275 * Si somma, trasposta, ad una matrice completa con metodi generici.
1276 */
AddToT(MatrixHandler & MH)1277 MatrixHandler& AddToT(MatrixHandler& MH) const {
1278 switch (eStatus) {
1279 case FULL:
1280 return FullSubMatrixHandler::AddToT(MH);
1281
1282 case SPARSE:
1283 return SparseSubMatrixHandler::AddToT(MH);
1284
1285 default:
1286 return MH;
1287 }
1288 };
1289
1290 /*
1291 * Si somma ad una matrice completa con metodi per matrici piene.
1292 */
AddTo(FullMatrixHandler & MH)1293 MatrixHandler& AddTo(FullMatrixHandler& MH) const {
1294 switch (eStatus) {
1295 case FULL:
1296 return FullSubMatrixHandler::AddTo(MH);
1297
1298 case SPARSE:
1299 return SparseSubMatrixHandler::AddTo(MH);
1300
1301 default:
1302 return MH;
1303 }
1304 };
1305
1306 /*
1307 * Si somma, trasposta, ad una matrice completa con metodi per matrici piene.
1308 */
AddToT(FullMatrixHandler & MH)1309 MatrixHandler& AddToT(FullMatrixHandler& MH) const {
1310 switch (eStatus) {
1311 case FULL:
1312 return FullSubMatrixHandler::AddToT(MH);
1313
1314 case SPARSE:
1315 return SparseSubMatrixHandler::AddToT(MH);
1316
1317 default:
1318 return MH;
1319 }
1320 };
1321
1322 /*
1323 * Si sottrae da una matrice completa con metodi generici.
1324 */
SubFrom(MatrixHandler & MH)1325 MatrixHandler& SubFrom(MatrixHandler& MH) const {
1326 switch (eStatus) {
1327 case FULL:
1328 return FullSubMatrixHandler::SubFrom(MH);
1329
1330 case SPARSE:
1331 return SparseSubMatrixHandler::SubFrom(MH);
1332
1333 default:
1334 return MH;
1335 }
1336 };
1337
1338 /*
1339 * Si sottrae, trasposta, da una matrice completa con metodi generici.
1340 */
SubFromT(MatrixHandler & MH)1341 MatrixHandler& SubFromT(MatrixHandler& MH) const {
1342 switch (eStatus) {
1343 case FULL:
1344 return FullSubMatrixHandler::SubFromT(MH);
1345
1346 case SPARSE:
1347 return SparseSubMatrixHandler::SubFromT(MH);
1348
1349 default:
1350 return MH;
1351 }
1352 };
1353
1354 /*
1355 * Si sottrae da una matrice completa con metodi per matrici piene.
1356 */
SubFrom(FullMatrixHandler & MH)1357 MatrixHandler& SubFrom(FullMatrixHandler& MH) const {
1358 switch (eStatus) {
1359 case FULL:
1360 return FullSubMatrixHandler::SubFrom(MH);
1361
1362 case SPARSE:
1363 return SparseSubMatrixHandler::SubFrom(MH);
1364
1365 default:
1366 return MH;
1367 }
1368 };
1369
1370 /*
1371 * Si sottrae, trasposta, da una matrice completa con metodi per matrici piene.
1372 */
SubFromT(FullMatrixHandler & MH)1373 MatrixHandler& SubFromT(FullMatrixHandler& MH) const {
1374 switch (eStatus) {
1375 case FULL:
1376 return FullSubMatrixHandler::SubFromT(MH);
1377
1378 case SPARSE:
1379 return SparseSubMatrixHandler::SubFromT(MH);
1380
1381 default:
1382 return MH;
1383 }
1384 };
1385
1386 const doublereal&
operator()1387 operator () (integer iRow, integer iCol) const {
1388 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1389 };
1390
1391 doublereal&
operator()1392 operator () (integer iRow, integer iCol) {
1393 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1394 };
1395 };
1396
1397 /* VariableSubMatrixHandler - end */
1398
1399
1400 /* SubVectorHandler - begin */
1401
1402 /*
1403 * Classe virtuale dei sottovettori.
1404 */
1405
1406 class SubVectorHandler : public VectorHandler {
1407 public:
1408 /* Costruttori */
1409
1410 /*
1411 * Distruttore virtuale.
1412 */
~SubVectorHandler(void)1413 virtual ~SubVectorHandler(void) {
1414 NO_OP;
1415 };
1416
1417 /* Metodi di servizio */
1418
1419 #ifdef DEBUG
1420 /*
1421 * Routine di verifica della validita' dell'oggetto.
1422 * Usata per il debug.
1423 */
1424 virtual void IsValid(void) const = 0;
1425 #endif /* DEBUG */
1426
1427 /* Operazioni su indici e coefficienti */
1428
1429 /*
1430 * Scrive un indice di riga
1431 */
1432 virtual void PutRowIndex(integer iSubRow, integer iRow) = 0;
1433
1434 /*
1435 * Ottiene un indice di riga
1436 */
1437 virtual integer iGetRowIndex(integer iSubRow) const = 0;
1438
1439 /*
1440 * Scrive una entry completa.
1441 * @param iSubRow numero della entry (indice del sotto-vettore)
1442 * @param iRow indice della entry
1443 * @param dCoef coefficiente della entry
1444 */
PutItem(integer iSubRow,integer iRow,const doublereal & dCoef)1445 virtual inline void PutItem(integer iSubRow, integer iRow,
1446 const doublereal& dCoef) {
1447 PutRowIndex(iSubRow, iRow);
1448 PutCoef(iSubRow, dCoef);
1449 };
1450
1451 /* Interazione con i vettori */
1452
1453 /*
1454 * Si somma ad un vettore con metodi generici
1455 */
1456 virtual VectorHandler& AddTo(VectorHandler& VH) const = 0;
1457 };
1458
1459 /*
1460 * Sottovettore standard, formato da un vettore di reali con associato
1461 * un vettore di interi che contiene gli indici di ogni coefficiente.
1462 * Per il vettore di reali viene usato un VectorHandler, da cui la classe e'
1463 * derivata.
1464 */
1465 class MySubVectorHandler : public SubVectorHandler, public MyVectorHandler {
1466 friend std::ostream&
1467 operator << (std::ostream& out, const SubVectorHandler& v);
1468
1469 protected:
1470 /* Puntatore all'array degli indici
1471 * Usato per rendere piu' efficiente l'accesso,
1472 * dato che gli indici sono a base 1, in stile FORTRAN
1473 */
1474 integer* piRowm1;
1475
1476 private:
1477 MySubVectorHandler(const MySubVectorHandler&);
1478
1479 public:
1480 /* Costruttori */
1481
1482 /*
1483 * Costruttore per memoria posseduta.
1484 * Specifica solo la dimensione dell'array, che deve essere non-nulla.
1485 * La memoria viene allocata e gestita dal VectorHandler.
1486 */
1487 MySubVectorHandler(integer iSize);
1488
1489 /*
1490 * Costruttore per memoria in prestito.
1491 * Riceve la dimensione dell'array e i puntatori alla memoria.
1492 */
1493 MySubVectorHandler(integer iSize, integer* piTmpRow,
1494 doublereal* pdTmpVec);
1495
1496 /*
1497 * Distruttore.
1498 */
~MySubVectorHandler(void)1499 virtual ~MySubVectorHandler(void) {
1500 Detach();
1501 };
1502
1503 /*
1504 * Tutti questi metodi sono richiesti perche'
1505 * la classe MySubVectorHandler dipende due volte
1506 * da VectorHandler e quindi vi e' un'ambiguita'
1507 * che va risolta (in realta' solo le funzioni di MyVectorHandler
1508 * sono definite, tutte le altre sono virtuali pure!)
1509 */
1510
1511 /* Metodi di servizio */
1512
1513 /*
1514 * Puntatore alla base del vettore (deprecato)
1515 */
pdGetVec(void)1516 virtual doublereal* pdGetVec(void) const {
1517 return MyVectorHandler::pdGetVec();
1518 };
1519
1520 /*
1521 * Dimensioni del vettore
1522 */
iGetSize(void)1523 virtual integer iGetSize(void) const {
1524 return MyVectorHandler::iGetSize();
1525 };
1526
1527 /*
1528 * Ridimensiona il vettore.
1529 * Nota: se il vettore possiede la memoria a cui punta,
1530 * la nuova dimensione puo' eccedere la massima dimensione corrente.
1531 */
1532 virtual void Resize(integer iSize);
1533
1534 /*
1535 * Inizializza il vettore con d
1536 */
Reset(void)1537 virtual void Reset(void) {
1538 MyVectorHandler::Reset();
1539 };
1540
1541 /*
1542 * Scollega il vettore dalla memoria ad esso associata.
1543 * Se il vettore possiede la memoria, viene delallocata.
1544 */
1545 void Detach(void);
1546
1547 /*
1548 * Collega il vettore alla memoria che gli viene passata.
1549 * La memoria a cui il vettore era collegato viene deallocata se
1550 * era posseduta dal vettore.
1551 */
1552 void
1553 Attach(integer iSize, doublereal* pd, integer* pi, integer iMSize = 0);
1554 #ifdef DEBUG
1555 /*
1556 * Verifica la validita' del vettore.
1557 * Usata per debug
1558 */
1559 virtual void IsValid(void) const;
1560 #endif /* DEBUG */
1561
1562 /* Operazioni sugli indici e sui coefficienti */
1563
1564 /*
1565 * Scrive un coefficiente in base al sottoindice.
1566 */
PutCoef(integer i,const doublereal & d)1567 virtual void PutCoef(integer i, const doublereal& d) {
1568 MyVectorHandler::PutCoef(i, d);
1569 };
1570
1571 /*
1572 * Incrementa un coefficiente in base al sottoindice.
1573 */
IncCoef(integer i,const doublereal & d)1574 virtual void IncCoef(integer i, const doublereal& d) {
1575 MyVectorHandler::IncCoef(i, d);
1576 };
1577
1578 /*
1579 * Decrementa un coefficiente in base al sottoindice.
1580 */
DecCoef(integer i,const doublereal & d)1581 virtual void DecCoef(integer i, const doublereal& d) {
1582 MyVectorHandler::DecCoef(i, d);
1583 };
1584
1585 /*
1586 * Ottiene un coefficiente in base al sottoindice.
1587 */
dGetCoef(integer i)1588 virtual const doublereal& dGetCoef(integer i) const {
1589 return MyVectorHandler::dGetCoef(i);
1590 };
1591
operator()1592 virtual inline const doublereal& operator () (integer iRow) const {
1593 return MyVectorHandler::operator () (iRow);
1594 };
1595
operator()1596 virtual inline doublereal& operator () (integer iRow) {
1597 return MyVectorHandler::operator () (iRow);
1598 };
1599
1600 /*
1601 * Scrive un indice di riga in base al sottoindice.
1602 */
1603 virtual inline void PutRowIndex(integer iSubRow, integer iRow);
1604
1605 /*
1606 * Ottiene un indice di riga in base al sottoindice.
1607 */
1608 virtual inline integer iGetRowIndex(integer iSubRow) const;
1609
1610 /*
1611 * Scrive una entry completa.
1612 * @param iSubRow numero della entry (indice del sotto-vettore)
1613 * @param iRow indice della entry
1614 * @param dCoef coefficiente della entry
1615 */
1616 virtual inline void
1617 PutItem(integer iSubRow, integer iRow, const doublereal& dCoef);
1618
1619 /* Interazione con i vettori */
1620
1621 /*
1622 * Si somma ad un vettore con metodi generici
1623 */
1624 virtual VectorHandler& AddTo(VectorHandler& VH) const;
1625
1626 /*
1627 * Si somma ad un MyVectorHandler
1628 */
1629 virtual VectorHandler& AddTo(MyVectorHandler& VH) const;
1630 };
1631
1632 inline void
PutRowIndex(integer iSubRow,integer iRow)1633 MySubVectorHandler::PutRowIndex(integer iSubRow, integer iRow)
1634 {
1635 #ifdef DEBUG
1636 IsValid();
1637 ASSERT((iSubRow > 0) && (iSubRow <= iCurSize));
1638 ASSERT(iRow > 0);
1639 #endif /* DEBUG */
1640
1641 piRowm1[iSubRow] = iRow;
1642 }
1643
1644 inline integer
iGetRowIndex(integer iSubRow)1645 MySubVectorHandler::iGetRowIndex(integer iSubRow) const
1646 {
1647 #ifdef DEBUG
1648 IsValid();
1649 ASSERT((iSubRow > 0) && (iSubRow <= iCurSize));
1650 #endif /* DEBUG */
1651
1652 return piRowm1[iSubRow];
1653 }
1654
1655 inline void
PutItem(integer iSubRow,integer iRow,const doublereal & dCoef)1656 MySubVectorHandler::PutItem(integer iSubRow, integer iRow,
1657 const doublereal& dCoef)
1658 {
1659 #ifdef DEBUG
1660 IsValid();
1661 ASSERT((iSubRow > 0) && (iSubRow <= iCurSize));
1662 ASSERT(iRow > 0);
1663 #endif /* DEBUG */
1664
1665 piRowm1[iSubRow] = iRow;
1666 pdVecm1[iSubRow] = dCoef;
1667 }
1668
1669 /* Operazioni esterne su SubMatrixHandler e su SubVectorHandler */
1670
1671 /*
1672 * Operatore per scrittura di SubVectorHandler su ostream.
1673 * Usato principalmente per debug
1674 */
1675 extern std::ostream&
1676 operator << (std::ostream& out, const SubVectorHandler& v);
1677
1678 /*
1679 * Operatore per scrittura di FullSubMatrixHandler su ostream.
1680 * Usato principalmente per debug
1681 */
1682 extern std::ostream&
1683 operator << (std::ostream& out, const FullSubMatrixHandler& m);
1684
1685
1686 /* SubVectorHandler - end */
1687
1688 #endif /* SUBMAT_H */
1689