MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
submat.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/submat.cc,v 1.49 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 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
35 
36 #include <string.h> /* for memset() */
37 #include <iomanip>
38 
39 #include <submat.h>
40 
41 /* SubMatrixHandler - begin */
42 
44 {
45  NO_OP;
46 }
47 
48 /* SubMatrixHandler - end */
49 
50 
51 /* FullSubMatrixHandler - begin */
52 
54  integer* piTmpVec,
55  integer iDoubleSize,
56  doublereal* pdTmpMat,
57  integer iMaxCols,
58  doublereal **ppdCols)
59 : FullMatrixHandler(pdTmpMat, ppdCols, iDoubleSize, 1, 1, iMaxCols),
60 iVecSize(iIntSize), piRowm1(0), piColm1(0)
61 {
62 #ifdef DEBUG
63  IsValid();
64 #endif
65 
66  /*
67  * Non posso fare controlli piu' precisi sulla consistenza
68  * delle dimensioni delle aree di lavoro perche' e' legale
69  * che la submatrix sia in seguito ridimensionata nei modi
70  * piu' vari. */
71  ASSERT(piTmpVec);
72  piRowm1 = piTmpVec - 1;
74 }
75 
77 : FullMatrixHandler(iNR, iNC), iVecSize(iNR + iNC), piRowm1(0), piColm1(0)
78 {
81 
82  piRowm1--;
83  piColm1 = piRowm1 + iNR;
84 }
85 
86 
88 {
89  if (bOwnsMemory) {
90  piRowm1++;
92  }
93 }
94 
95 #ifdef DEBUG
96 void
97 FullSubMatrixHandler::IsValid(void) const
98 {
99  FullMatrixHandler::IsValid();
100 
101  ASSERT(iVecSize > 0);
103  ASSERT(piRowm1 != NULL);
104  ASSERT(piColm1 != NULL);
105 
106 #ifdef DEBUG_MEMMANAGER
107  ASSERT(defaultMemoryManager.fIsValid(piRowm1 + 1,
108  iVecSize*sizeof(integer)));
109 #endif /* DEBUG_MEMMANAGER */
110 }
111 #endif /* DEBUG */
112 
113 
114 void
116 {
117 #ifdef DEBUG
118  IsValid();
119 #endif /* DEBUG */
120 
122 
123 #if 0
124  /*
125  * this is not strictly required, because all the indices should
126  * be explicitly set before the matrix is used
127  */
128  for (integer i = iNumRows + iNumCols; i > 0; i--) {
129  piRowm1[i] = 0;
130  }
131 #endif
132 }
133 
134 /*
135  * Modifica le dimensioni correnti
136  */
137 void
139 #ifdef DEBUG
140  IsValid();
141 #endif /* DEBUG */
142 
143  ASSERT(iNewRow > 0);
144  ASSERT(iNewCol > 0);
145  ASSERT(iNewRow + iNewCol <= iVecSize);
146  ASSERT(iNewRow*iNewCol <= iRawSize);
147 
148  if (iNewRow <= 0
149  || iNewCol <= 0
150  || iNewRow + iNewCol > iVecSize
151  || iNewRow*iNewCol > iRawSize
152  || iNewCol > iMaxCols) {
153  silent_cerr("FullSubMatrixHandler::Resize() - error"
154  << std::endl);
155 
157  }
158 
159  iNumRows = iNewRow;
160  iNumCols = iNewCol;
161 
162  /*
163  * FIXME: don't call FullMatrixHandler::Resize, because
164  * we know there's room enough and we don't want to
165  * waste time regenerating support stuff; we simply
166  * use a portion of the matrix as is.
167  */
168 #if 0
169  FullMatrixHandler::Resize(iNewRow, iNewCol);
170 
171  ASSERT(piRowm1 != NULL);
172  piColm1 = piRowm1 + iNewRow;
173 #endif
174 
175 #ifdef DEBUG
176  IsValid();
177 #endif /* DEBUG */
178 }
179 
180 /* Ridimensiona ed inizializza. */
181 void
183 {
186 }
187 
188 /*
189  * Collega la matrice Full alla memoria che gli viene passata
190  * in ingresso
191  */
192 void
193 FullSubMatrixHandler::Attach(int iRows, int iCols, integer* piTmpIndx)
194 {
195  iNumRows = iRows;
196  iNumCols = iCols;
197 
198  if (bOwnsMemory) {
199  if (piRowm1 != 0) {
200  piRowm1++;
202  }
203  piRowm1 = piTmpIndx - 1;
205  }
206 }
207 
208 void
210 {
211  FullMatrixHandler::Add(iRow, iCol, v);
212 }
213 
214 void
216 {
217  FullMatrixHandler::Sub(iRow, iCol, v);
218 }
219 
220 void
222 {
223  FullMatrixHandler::Put(iRow, iCol, v);
224 }
225 
226 void
228 {
229  FullMatrixHandler::AddT(iRow, iCol, v);
230 }
231 
232 void
234 {
235  FullMatrixHandler::SubT(iRow, iCol, v);
236 }
237 
238 void
240 {
241  FullMatrixHandler::PutT(iRow, iCol, v);
242 }
243 
244 #if 0 /* FIXME: replace original? */
245 /* somma un vettore di tipo Vec3 in una data posizione in diagonale */
246 void
247 FullSubMatrixHandler::AddDiag(integer iRow, integer iCol, const Vec3& v)
248 {
249  /* iRow e iCol sono gli indici effettivi di riga e colonna
250  * es. per il primo coefficiente:
251  * iRow = 1, iCol = 1 */
252 
253 #ifdef DEBUG
254  IsValid();
255 
256  ASSERT(iRow > 0);
257  ASSERT(iRow <= iNumRows - 2);
258  ASSERT(iCol > 0);
259  ASSERT(iCol <= iNumCols - 2);
260 #endif /* DEBUG */
261 
262  /* Nota: assume che la matrice sia organizzata
263  * "per colonne" (stile FORTRAN)
264  */
265  doublereal* pdFrom = v.pGetVec();
266 
267  ppdColsm1[iCol][iRow] += pdFrom[V1];
268 
269  ppdColsm1[iCol + 1][iRow + 1] += pdFrom[V2];
270 
271  ppdColsm1[iCol + 2][iRow + 2] += pdFrom[V3];
272 }
273 
274 /* sottrae un vettore di tipo Vec3 in una data posizione in diagonale */
275 void
276 FullSubMatrixHandler::SubDiag(integer iRow, integer iCol, const Vec3& v)
277 {
278  /* iRow e iCol sono gli indici effettivi di riga e colonna
279  * es. per il primo coefficiente:
280  * iRow = 1, iCol = 1 */
281 
282 #ifdef DEBUG
283  IsValid();
284 
285  ASSERT(iRow > 0);
286  ASSERT(iRow <= iNumRows - 2);
287  ASSERT(iCol > 0);
288  ASSERT(iCol <= iNumCols - 2);
289 #endif /* DEBUG */
290 
291  /* Nota: assume che la matrice sia organizzata
292  * "per colonne" (stile FORTRAN)
293  */
294  doublereal* pdFrom = v.pGetVec();
295 
296  ppdColsm1[iCol][iRow] -= pdFrom[V1];
297 
298  ppdColsm1[iCol + 1][iRow + 1] -= pdFrom[V2];
299 
300  ppdColsm1[iCol + 2][iRow + 2] -= pdFrom[V3];
301 }
302 
303 /* scrive un vettore di tipo Vec3 in una data posizione in diagonale */
304 void
306 {
307  /* iRow e iCol sono gli indici effettivi di riga e colonna
308  * es. per il primo coefficiente:
309  * iRow = 1, iCol = 1 */
310 
311 #ifdef DEBUG
312  IsValid();
313 
314  ASSERT(iRow > 0);
315  ASSERT(iRow <= iNumRows - 2);
316  ASSERT(iCol > 0);
317  ASSERT(iCol <= iNumCols - 2);
318 #endif /* DEBUG */
319 
320  /* Nota: assume che la matrice sia organizzata
321  * "per colonne" (stile FORTRAN)
322  */
323  doublereal* pdFrom = v.pGetVec();
324 
325  ppdColsm1[iCol][iRow] = pdFrom[V1];
326 
327  ppdColsm1[iCol + 1][iRow + 1] = pdFrom[V2];
328 
329  ppdColsm1[iCol + 2][iRow + 2] = pdFrom[V3];
330 }
331 
332 /* somma un vettore di tipo Vec3 in una data posizione [ v x ] */
333 void
334 FullSubMatrixHandler::AddCross(integer iRow, integer iCol, const Vec3& v)
335 {
336  /* iRow e iCol sono gli indici effettivi di riga e colonna
337  * es. per il primo coefficiente:
338  * iRow = 1, iCol = 1 */
339 
340 #ifdef DEBUG
341  IsValid();
342 
343  ASSERT(iRow > 0);
344  ASSERT(iRow <= iNumRows - 2);
345  ASSERT(iCol > 0);
346  ASSERT(iCol <= iNumCols - 2);
347 #endif /* DEBUG */
348 
349  /* Nota: assume che la matrice sia organizzata
350  * "per colonne" (stile FORTRAN)
351  */
352  doublereal* pdFrom = v.pGetVec();
353 
354  ppdColsm1[iCol + 1][iRow] -= pdFrom[V3];
355  ppdColsm1[iCol + 2][iRow] += pdFrom[V2];
356 
357  ppdColsm1[iCol][iRow + 1] += pdFrom[V3];
358  ppdColsm1[iCol + 2][iRow + 1] -= pdFrom[V1];
359 
360  ppdColsm1[iCol][iRow + 2] -= pdFrom[V2];
361  ppdColsm1[iCol + 1][iRow + 2] += pdFrom[V1];
362 }
363 
364 /* sottrae un vettore di tipo Vec3 in una data posizione [ v x ] */
365 void
366 FullSubMatrixHandler::SubCross(integer iRow, integer iCol, const Vec3& v)
367 {
368  /* iRow e iCol sono gli indici effettivi di riga e colonna
369  * es. per il primo coefficiente:
370  * iRow = 1, iCol = 1 */
371 
372 #ifdef DEBUG
373  IsValid();
374 
375  ASSERT(iRow > 0);
376  ASSERT(iRow <= iNumRows - 2);
377  ASSERT(iCol > 0);
378  ASSERT(iCol <= iNumCols - 2);
379 #endif /* DEBUG */
380 
381  /* Nota: assume che la matrice sia organizzata
382  * "per colonne" (stile FORTRAN)
383  */
384  doublereal* pdFrom = v.pGetVec();
385 
386  ppdColsm1[iCol + 1][iRow] += pdFrom[V3];
387  ppdColsm1[iCol + 2][iRow] -= pdFrom[V2];
388 
389  ppdColsm1[iCol][iRow + 1] -= pdFrom[V3];
390  ppdColsm1[iCol + 2][iRow + 1] += pdFrom[V1];
391 
392  ppdColsm1[iCol][iRow + 2] += pdFrom[V2];
393  ppdColsm1[iCol + 1][iRow + 2] -= pdFrom[V1];
394 }
395 
396 /* scrive un vettore di tipo Vec3 in una data posizione [ v x ] */
397 void
399 {
400  /* iRow e iCol sono gli indici effettivi di riga e colonna
401  * es. per il primo coefficiente:
402  * iRow = 1, iCol = 1 */
403 
404 #ifdef DEBUG
405  IsValid();
406 
407  ASSERT(iRow > 0);
408  ASSERT(iRow <= iNumRows - 2);
409  ASSERT(iCol > 0);
410  ASSERT(iCol <= iNumCols - 2);
411 #endif /* DEBUG */
412 
413  /* Nota: assume che la matrice sia organizzata
414  * "per colonne" (stile FORTRAN)
415  */
416  doublereal* pdFrom = v.pGetVec();
417 
418  ppdColsm1[iCol + 1][iRow] = -pdFrom[V3];
419  ppdColsm1[iCol + 2][iRow] = pdFrom[V2];
420 
421  ppdColsm1[iCol][iRow + 1] = pdFrom[V3];
422  ppdColsm1[iCol + 2][iRow + 1] = -pdFrom[V1];
423 
424  ppdColsm1[iCol][iRow + 2] = -pdFrom[V2];
425  ppdColsm1[iCol + 1][iRow + 2] = pdFrom[V1];
426 }
427 #endif
428 
429 void
431 {
432  FullMatrixHandler::Add(iRow, iCol, m);
433 }
434 
435 void
437 {
438  FullMatrixHandler::AddT(iRow, iCol, m);
439 }
440 
441 void
443 {
444  FullMatrixHandler::Sub(iRow, iCol, m);
445 }
446 
447 void
449 {
450  FullMatrixHandler::SubT(iRow, iCol, m);
451 }
452 
453 void
455 {
456  FullMatrixHandler::Put(iRow, iCol, m);
457 }
458 
459 void
461 {
462  FullMatrixHandler::PutT(iRow, iCol, m);
463 }
464 
465 void
467 {
468  FullMatrixHandler::Add(iRow, iCol, m);
469 }
470 
471 void
473 {
474  FullMatrixHandler::Sub(iRow, iCol, m);
475 }
476 
477 void
479 {
480  FullMatrixHandler::Put(iRow, iCol, m);
481 }
482 
483 void
485 {
486  FullMatrixHandler::AddT(iRow, iCol, m);
487 }
488 
489 void
491 {
492  FullMatrixHandler::SubT(iRow, iCol, m);
493 }
494 
495 void
497 {
498  FullMatrixHandler::PutT(iRow, iCol, m);
499 }
500 
501 void
503 {
504  FullMatrixHandler::Add(iRow, iCol, m);
505 }
506 
507 /* sottrae una matrice di tipo MatNx3 in una data posizione */
508 void
510 {
511  FullMatrixHandler::Sub(iRow, iCol, m);
512 }
513 
514 /* setta una matrice di tipo MatNx3 in una data posizione */
515 void
517 {
518  FullMatrixHandler::Put(iRow, iCol, m);
519 }
520 
521 /* setta una matrice di tipo Mat3xN in una data posizione */
522 void
524  const Vec3& v)
525 {
526  /* iFirstRow e iFirstCol sono gli indici effettivi di riga e colonna -1
527  * es. per il primo coefficiente:
528  * iRow = 0, iCol = 0 */
529 
530 #ifdef DEBUG
531  IsValid();
532 
533  ASSERT(iFirstRow >= 0);
534  ASSERT(iFirstRow <= iNumRows - 3);
535  ASSERT(iFirstCol >= 0);
536  ASSERT(iFirstCol <= iNumCols - 3);
537 #endif /* DEBUG */
538 
539  const doublereal *pdv = v.pGetVec();
540 
541  ppdColsm1[iFirstCol + 1][iFirstRow + 1] = pdv[V1];
542  ppdColsm1[iFirstCol + 2][iFirstRow + 2] = pdv[V2];
543  ppdColsm1[iFirstCol + 3][iFirstRow + 3] = pdv[V3];
544 }
545 
546 
547 /* setta una matrice di tipo Mat3xN in una data posizione */
548 void
550  const doublereal& d)
551 {
552  /* iFirstRow e iFirstCol sono gli indici effettivi di riga e colonna -1
553  * es. per il primo coefficiente:
554  * iRow = 0, iCol = 0 */
555 
556 #ifdef DEBUG
557  IsValid();
558 
559  ASSERT(iFirstRow >= 0);
560  ASSERT(iFirstRow <= iNumRows - 3);
561  ASSERT(iFirstCol >= 0);
562  ASSERT(iFirstCol <= iNumCols - 3);
563 #endif /* DEBUG */
564 
565  ppdColsm1[iFirstCol + 1][iFirstRow + 1] = d;
566  ppdColsm1[iFirstCol + 2][iFirstRow + 2] = d;
567  ppdColsm1[iFirstCol + 3][iFirstRow + 3] = d;
568 }
569 
570 
571 /* setta una matrice di tipo Mat3xN in una data posizione */
572 void
574  const Vec3& v)
575 {
576  /* iFirstRow e iFirstCol sono gli indici effettivi di riga e colonna -1
577  * es. per il primo coefficiente:
578  * iRow = 0, iCol = 0 */
579 
580 #ifdef DEBUG
581  IsValid();
582 
583  ASSERT(iFirstRow >= 0);
584  ASSERT(iFirstRow <= iNumRows - 3);
585  ASSERT(iFirstCol >= 0);
586  ASSERT(iFirstCol <= iNumCols - 3);
587 #endif /* DEBUG */
588 
589  const doublereal *pdv = v.pGetVec();
590 
591  ppdColsm1[iFirstCol + 1][ iFirstRow + 2] = -pdv[V3];
592  ppdColsm1[iFirstCol + 1][ iFirstRow + 3] = pdv[V2];
593 
594  ppdColsm1[iFirstCol + 2][ iFirstRow + 1] = pdv[V3];
595  ppdColsm1[iFirstCol + 2][ iFirstRow + 3] = -pdv[V1];
596 
597  ppdColsm1[iFirstCol + 3][ iFirstRow + 1] = -pdv[V2];
598  ppdColsm1[iFirstCol + 3][ iFirstRow + 2] = pdv[V1];
599 }
600 
601 
602 /* somma la matrice ad un matrix handler usando i metodi generici */
605 {
606 #ifdef DEBUG
607  IsValid();
608  MH.IsValid();
609 #endif /* DEBUG */
610 
611  ASSERT(MH.iGetNumRows() >= iNumRows);
612  ASSERT(MH.iGetNumCols() >= iNumCols);
613 
614  for (integer c = iNumCols; c > 0; c--) {
615  ASSERT(piColm1[c] > 0);
616  ASSERT(piColm1[c] <= MH.iGetNumCols());
617 
618  for (integer r = iNumRows; r > 0; r--) {
619  ASSERT(piRowm1[r] > 0);
620  ASSERT(piRowm1[r] <= MH.iGetNumRows());
621 
622  MH(piRowm1[r], piColm1[c]) += ppdColsm1[c][r];
623  }
624  }
625 
626  return MH;
627 }
628 
629 
630 /* somma la matrice, trasposta, ad un matrix handler usando i metodi generici */
633 {
634 #ifdef DEBUG
635  IsValid();
636  MH.IsValid();
637 #endif /* DEBUG */
638 
639  ASSERT(MH.iGetNumRows() >= iNumCols);
640  ASSERT(MH.iGetNumCols() >= iNumRows);
641 
642  for (integer c = iNumCols; c > 0; c--) {
643  ASSERT(piColm1[c] > 0);
644  ASSERT(piColm1[c] <= MH.iGetNumRows());
645 
646  for (integer r = iNumRows; r > 0; r--) {
647  ASSERT(piRowm1[r] > 0);
648  ASSERT(piRowm1[r] <= MH.iGetNumCols());
649 
650  MH(piColm1[c], piRowm1[r]) += ppdColsm1[c][r];
651  }
652  }
653 
654  return MH;
655 }
656 
657 
658 /* somma la matrice ad un FullMatrixHandler */
661 {
662 #ifdef DEBUG
663  IsValid();
664  MH.IsValid();
665 #endif /* DEBUG */
666 
667  ASSERT(MH.iGetNumRows() >= iNumRows);
668  ASSERT(MH.iGetNumCols() >= iNumCols);
669 
670  doublereal **ppd = MH.ppdColsm1;
671 
672  for (integer c = iNumCols; c > 0; c--) {
673  ASSERT(piColm1[c] > 0);
674  ASSERT(piColm1[c] <= MH.iGetNumCols());
675 
676  for (integer r = iNumRows; r > 0; r--) {
677  ASSERT(piRowm1[r] > 0);
678  ASSERT(piRowm1[r] <= MH.iGetNumRows());
679 
680  ppd[piColm1[c]][piRowm1[r]] += ppdColsm1[c][r];
681  }
682  }
683 
684  return MH;
685 }
686 
687 
688 /* somma la matrice, trasposta, ad un FullMatrixHandler */
691 {
692 #ifdef DEBUG
693  IsValid();
694  MH.IsValid();
695 #endif /* DEBUG */
696 
697  ASSERT(MH.iGetNumRows() >= iNumCols);
698  ASSERT(MH.iGetNumCols() >= iNumRows);
699 
700  doublereal **ppd = MH.ppdColsm1;
701 
702  for (integer c = iNumCols; c > 0; c--) {
703  ASSERT(piColm1[c] > 0);
704  ASSERT(piColm1[c] <= MH.iGetNumRows());
705 
706  for (integer r = iNumRows; r > 0; r--) {
707  ASSERT(piRowm1[r] > 0);
708  ASSERT(piRowm1[r] <= MH.iGetNumCols());
709 
710  ppd[piRowm1[r]][piColm1[c]] += ppdColsm1[c][r];
711  }
712  }
713 
714  return MH;
715 }
716 
717 
718 /* sottrae la matrice da un matrix handler usando i metodi generici */
721 {
722 #ifdef DEBUG
723  IsValid();
724  MH.IsValid();
725 #endif /* DEBUG */
726 
727  ASSERT(MH.iGetNumRows() >= iNumRows);
728  ASSERT(MH.iGetNumCols() >= iNumCols);
729 
730  for (integer c = iNumCols; c > 0; c--) {
731  ASSERT(piColm1[c] > 0);
732  ASSERT(piColm1[c] <= MH.iGetNumCols());
733 
734  for (integer r = iNumRows; r > 0; r--) {
735  ASSERT(piRowm1[r] > 0);
736  ASSERT(piRowm1[r] <= MH.iGetNumRows());
737 
738  MH(piRowm1[r], piColm1[c]) -= ppdColsm1[c][r];
739  }
740  }
741 
742  return MH;
743 }
744 
745 
746 /* sottrae la matrice, trasposta, da un matrix handler usando i metodi generici */
749 {
750 #ifdef DEBUG
751  IsValid();
752  MH.IsValid();
753 #endif /* DEBUG */
754 
755  ASSERT(MH.iGetNumRows() >= iNumCols);
756  ASSERT(MH.iGetNumCols() >= iNumRows);
757 
758  for (integer c = iNumCols; c > 0; c--) {
759  ASSERT(piColm1[c] > 0);
760  ASSERT(piColm1[c] <= MH.iGetNumRows());
761 
762  for (integer r = iNumRows; r > 0; r--) {
763  ASSERT(piRowm1[r] > 0);
764  ASSERT(piRowm1[r] <= MH.iGetNumCols());
765 
766  MH(piColm1[c], piRowm1[r]) -= ppdColsm1[c][r];
767  }
768  }
769 
770  return MH;
771 }
772 
773 
774 /* sottrae la matrice da un FullMatrixHandler */
777 {
778 #ifdef DEBUG
779  IsValid();
780  MH.IsValid();
781 #endif /* DEBUG */
782 
783  ASSERT(MH.iGetNumRows() >= iNumRows);
784  ASSERT(MH.iGetNumCols() >= iNumCols);
785 
786  doublereal **ppd = MH.ppdColsm1;
787 
788  for (integer c = iNumCols; c > 0; c--) {
789  ASSERT(piColm1[c] > 0);
790  ASSERT(piColm1[c] <= MH.iGetNumCols());
791 
792  for (integer r = iNumRows; r > 0; r--) {
793  ASSERT(piRowm1[r] > 0);
794  ASSERT(piRowm1[r] <= MH.iGetNumRows());
795 
796  ppd[piColm1[c]][piRowm1[r]] -= ppdColsm1[c][r];
797  }
798  }
799 
800  return MH;
801 };
802 
803 
804 /* sottrae la matrice, trasposta, da un FullMatrixHandler */
807 {
808 #ifdef DEBUG
809  IsValid();
810  MH.IsValid();
811 #endif /* DEBUG */
812 
813  ASSERT(MH.iGetNumRows() >= iNumCols);
814  ASSERT(MH.iGetNumCols() >= iNumRows);
815 
816  doublereal **ppd = MH.ppdColsm1;
817 
818  for (integer c = iNumCols; c > 0; c--) {
819  ASSERT(piColm1[c] > 0);
820  ASSERT(piColm1[c] <= MH.iGetNumRows());
821 
822  for (integer r = iNumRows; r > 0; r--) {
823  ASSERT(piRowm1[r] > 0);
824  ASSERT(piRowm1[r] <= MH.iGetNumCols());
825 
826  ppd[piRowm1[r]][piColm1[c]] -= ppdColsm1[c][r];
827  }
828  }
829 
830  return MH;
831 };
832 
833 
834 /* output, usato principalmente per debug */
835 std::ostream&
836 operator << (std::ostream& out, const FullSubMatrixHandler& m)
837 {
838 #ifdef DEBUG
839  m.IsValid();
840 #endif /* DEBUG */
841 
842  ASSERT(m.iNumRows > 0);
843  ASSERT(m.iNumCols > 0);
844  ASSERT(m.piRowm1 != NULL);
845  ASSERT(m.piColm1 != NULL);
846  ASSERT(m.ppdColsm1 != NULL);
847 
848  out << std::setw(12) << "";
849  for (integer c = 1; c <= m.iNumCols; c++) {
850  out << std::setw(12) << m.piColm1[c];
851  }
852  out << std::endl << std::endl;
853 
854  for (integer r = 1; r <= m.iNumRows; r++) {
855  out << std::setw(12) << m.piRowm1[r];
856  for (integer c = 1; c <= m.iNumCols; c++) {
857  out << std::setw(12) << m.ppdColsm1[c][r];
858  }
859  out << std::endl;
860  }
861 
862  return out << std::endl;
863 }
864 
865 void
867  integer iCol,
868  const FullMatrixHandler & source)
869 {
870  FullMatrixHandler::Add(iRow, iCol, source);
871 }
872 
873 void
875  integer iCol,
876  const FullMatrixHandler & source)
877 {
878  FullMatrixHandler::Sub(iRow, iCol, source);
879 }
880 
881 void
883  integer iCol,
884  const FullMatrixHandler & source)
885 {
886  FullMatrixHandler::Put(iRow, iCol, source);
887 }
888 
889 void
891  integer iCol,
892  const FullMatrixHandler & source,
893  const doublereal dCoef)
894 {
895  FullMatrixHandler::Add(iRow, iCol, source, dCoef);
896 }
897 
898 void
900  integer iCol,
901  const FullMatrixHandler & source,
902  const doublereal dCoef)
903 {
904  FullMatrixHandler::Sub(iRow, iCol, source, dCoef);
905 }
906 
907 void
909  integer iCol,
910  const FullMatrixHandler & source,
911  const doublereal dCoef)
912 {
913  FullMatrixHandler::Put(iRow, iCol, source, dCoef);
914 }
915 
916 void
918  integer iCol,
919  const FullMatrixHandler & source)
920 {
921  FullMatrixHandler::AddT(iRow, iCol, source);
922 }
923 
924 void
926  integer iCol,
927  const FullMatrixHandler & source)
928 {
929  FullMatrixHandler::SubT(iRow, iCol, source);
930 }
931 
932 void
934  integer iCol,
935  const FullMatrixHandler & source)
936 {
937  FullMatrixHandler::PutT(iRow, iCol, source);
938 }
939 
940 void
942  integer iCol,
943  const FullMatrixHandler & source,
944  const doublereal dCoef)
945 {
946  FullMatrixHandler::AddT(iRow, iCol, source, dCoef);
947 }
948 
949 void
951  integer iCol,
952  const FullMatrixHandler & source,
953  const doublereal dCoef)
954 {
955  FullMatrixHandler::SubT(iRow, iCol, source, dCoef);
956 }
957 
958 void
960  integer iCol,
961  const FullMatrixHandler & source,
962  const doublereal dCoef)
963 {
964  FullMatrixHandler::PutT(iRow, iCol, source, dCoef);
965 }
966 
967 /* FullSubMatrixHandler - end */
968 
969 
970 /* SparseSubMatrixHandler - begin */
971 
973  integer* piTmpIndex, integer iTmpDouble, doublereal* pdTmpMat)
974 : bOwnsMemory(false),
975 iIntSize(iTmpInt), iDoubleSize(iTmpDouble),
976 iNumItems(iTmpInt/2),
977 piRowm1(0), piColm1(0),
978 pdMatm1(0)
979 {
980  ASSERT(piTmpIndex);
981  ASSERT(pdTmpMat);
982 
983  piRowm1 = piTmpIndex - 1;
985 
986  pdMatm1 = pdTmpMat - 1;
987 
988 #ifdef DEBUG
989  IsValid();
990 #endif /* DEBUG */
991 }
992 
994 : iIntSize(2*iTmpInt), iDoubleSize(iTmpInt),
995 iNumItems(iTmpInt), piRowm1(0), piColm1(0), pdMatm1(0)
996 {
998  pdMatm1--;
999 
1001  piRowm1--;
1003 
1004  bOwnsMemory = true;
1005 }
1006 
1007 /* Distruttore banale.
1008  * Nota: dato che la classe non possiede la memoria,
1009  * non ne deve deallocare
1010  */
1012 {
1013  if (bOwnsMemory) {
1014  pdMatm1++;
1016 
1017  piRowm1++;
1019  }
1020 }
1021 
1022 #ifdef DEBUG
1023 void
1024 SparseSubMatrixHandler::IsValid(void) const
1025 {
1026  ASSERT(iIntSize > 0);
1027  ASSERT(iIntSize%2 == 0);
1028  ASSERT(iDoubleSize > 0);
1029  ASSERT(iIntSize >= 2*iDoubleSize);
1030  ASSERT(iNumItems >= 0);
1031  ASSERT(piRowm1 != NULL);
1032  ASSERT(piColm1 != NULL);
1033  ASSERT(pdMatm1 != NULL);
1034 
1035 #ifdef DEBUG_MEMMANAGER
1036  ASSERT(defaultMemoryManager.fIsValid(piRowm1 + 1,
1037  iIntSize*sizeof(integer)));
1038  ASSERT(defaultMemoryManager.fIsValid(pdMatm1 + 1,
1039  iDoubleSize*sizeof(doublereal)));
1040 #endif /* DEBUG_MEMMANAGER */
1041 }
1042 #endif /* DEBUG */
1043 
1044 
1045 /*
1046  * Ridimensiona la matrice.
1047  * Nota: solo il primo argomento viene considerato,
1048  * e rappresenta il numero totale di entries.
1049  * Questo metodo deve essere chiamato prima di qualsiasi
1050  * operazione sulla matrice.
1051  */
1052 void
1054 {
1055 #ifdef DEBUG
1056  IsValid();
1057 #endif /* DEBUG */
1058 
1059  ASSERT(iNewRow > 0);
1060  ASSERT(2*iNewRow <= iIntSize);
1061  ASSERT(iNewRow <= iDoubleSize);
1062  ASSERT(piRowm1 != NULL);
1063 
1064  if (iNewRow <= 0
1065  || 2*iNewRow > iIntSize
1066  || iNewRow > iDoubleSize) {
1067  silent_cerr("SparseSubMatrixHandler::Resize() - error"
1068  << std::endl);
1070  }
1071 
1072  iNumItems = iNewRow;
1073 
1074 #ifdef DEBUG
1075  IsValid();
1076 #endif /* DEBUG */
1077 }
1078 
1079 /*
1080  * Ridimensiona ed inizializza.
1081  * Unione dei due metodi precedenti
1082  */
1083 void
1085 {
1086  Resize(iNewRow, iNewCol);
1087  Reset();
1088 }
1089 
1090 /*
1091  * Collega la matrice sparsa alla memoria che gli viene passata
1092  * in ingresso
1093  */
1094 void
1096  integer* piTmpIndx)
1097 {
1098  if (bOwnsMemory) {
1099  piRowm1++;
1101 
1102  pdMatm1++;
1104 
1105  bOwnsMemory = false;
1106  }
1107 
1108  ASSERT(iNumEntr > 0);
1109  ASSERT(piTmpIndx);
1110  ASSERT(pdTmpMat);
1111 
1112  iIntSize = iNumEntr*2;
1113  iDoubleSize = iNumEntr;
1114  iNumItems = iNumEntr;
1115  piRowm1 = piTmpIndx - 1;
1116  piColm1 = piRowm1 + iNumEntr;
1117  pdMatm1 = pdTmpMat - 1;
1118 
1119 #ifdef DEBUG
1120  IsValid();
1121 #endif /* DEBUG */
1122 }
1123 
1124 void
1126  integer iFirstCol, const Vec3& v)
1127 {
1128 #ifdef DEBUG
1129  IsValid();
1130 
1131  ASSERT(iNumItems >= 3);
1132  ASSERT(iSubIt > 0);
1133  ASSERT(iSubIt <= iNumItems - 2);
1134  ASSERT(iFirstRow >= 0);
1135  ASSERT(iFirstCol >= 0);
1136 #endif /* DEBUG */
1137 
1138  /*
1139  * Attenzione agli argomenti:
1140  * iSubIt e' il primo indice della matrice da utilizzare,
1141  * con 1 <= iSubit <= iCurSize;
1142  * iFirstRow e' il primo indice di riga -1, ovvero il primo indice
1143  * di riga della sottomatrice diag(v) piena e' iFirstRow + 1
1144  * iFirstCol e' il primo indice di colonna -1, ovvero il primo indice
1145  * di colonna della sottomatrice diag(v) piena e' iFirstCol + 1
1146  * v e' il vettore che genera diag(v)
1147  */
1148 
1149  /* Matrice diag(v) :
1150  *
1151  * 1 2 3
1152  *
1153  * 1 | v1 0 0 |
1154  * 2 | 0 v2 0 |
1155  * 3 | 0 0 v3 |
1156  */
1157 
1158  /* assume che il Vec3 sia un'array di 3 reali */
1159  const doublereal* pdFrom = v.pGetVec();
1160 
1161  doublereal* pdm = pdMatm1 + iSubIt;
1162  integer* pir = piRowm1 + iSubIt;
1163  integer* pic = piColm1 + iSubIt;
1164 
1165  /* Coefficiente 1,1 */
1166  pdm[0] = pdFrom[V1];
1167  pir[0] = iFirstRow + 1;
1168  pic[0] = iFirstCol + 1;
1169 
1170  /* Coefficiente 2,2 */
1171  pdm[1] = pdFrom[V2];
1172  pir[1] = iFirstRow + 2;
1173  pic[1] = iFirstCol + 2;
1174 
1175  /* Coefficiente 3,3 */
1176  pdm[2] = pdFrom[V3];
1177  pir[2] = iFirstRow + 3;
1178  pic[2] = iFirstCol + 3;
1179 }
1180 
1181 
1182 void
1184  integer iFirstCol, const doublereal& d)
1185 {
1186 #ifdef DEBUG
1187  IsValid();
1188 
1189  ASSERT(iNumItems >= 3);
1190  ASSERT(iSubIt > 0);
1191  ASSERT(iSubIt <= iNumItems - 2);
1192  ASSERT(iFirstRow >= 0);
1193  ASSERT(iFirstCol >= 0);
1194 #endif /* DEBUG */
1195 
1196  /* Attenzione agli argomenti:
1197  * iSubIt e' il primo indice della matrice da utilizzare,
1198  * con 1 <= iSubit <= iCurSize;
1199  * iFirstRow e' il primo indice di riga -1, ovvero il
1200  * primo indice di riga della sottomatrice I*d piena e' iFirstRow+1
1201  * iFirstCol e' il primo indice di colonna -1, ovvero il
1202  * primo indice di colonna della sottomatrice I*d piena e' iFirstCol+1
1203  * v e' il vettore che genera I*d */
1204 
1205  /* Matrice I*d :
1206  *
1207  * 1 2 3
1208  *
1209  * 1 | d 0 0 |
1210  * 2 | 0 d 0 |
1211  * 3 | 0 0 d |
1212  */
1213 
1214  doublereal* pdm = pdMatm1 + iSubIt;
1215  integer* pir = piRowm1 + iSubIt;
1216  integer* pic = piColm1 + iSubIt;
1217 
1218  /* Coefficiente 1,1 */
1219  pdm[0] = d;
1220  pir[0] = iFirstRow + 1;
1221  pic[0] = iFirstCol + 1;
1222 
1223  /* Coefficiente 2,2 */
1224  pdm[1] = d;
1225  pir[1] = iFirstRow + 2;
1226  pic[1] = iFirstCol + 2;
1227 
1228  /* Coefficiente 3,3 */
1229  pdm[2] = d;
1230  pir[2] = iFirstRow + 3;
1231  pic[2] = iFirstCol + 3;
1232 }
1233 
1234 
1235 void
1237  integer iFirstCol, const Vec3& v)
1238 {
1239 #ifdef DEBUG
1240  IsValid();
1241 
1242  ASSERT(iNumItems >= 6);
1243  ASSERT(iSubIt > 0);
1244  ASSERT(iSubIt <= iNumItems - 5);
1245  ASSERT(iFirstRow >= 0);
1246  ASSERT(iFirstCol >= 0);
1247 #endif /* DEBUG */
1248 
1249  /* Attenzione agli argomenti:
1250  * iSubIt e' il primo indice della matrice da utilizzare,
1251  * con 1 <= iSubit <= iCurSize;
1252  * iFirstRow e' il primo indice di riga -1, ovvero il
1253  * primo indice di riga della sottomatrice v/\ piena e' iFirstRow+1
1254  * iFirstCol e' il primo indice di colonna -1, ovvero il
1255  * primo indice di colonna della sottomatrice v/\ piena e' iFirstCol+1
1256  * v e' il vettore che genera v/\ */
1257 
1258  /* Matrice v/\ :
1259  *
1260  * 1 2 3
1261  *
1262  * 1 | 0 -v3 v2 |
1263  * 2 | v3 0 -v1 |
1264  * 3 | -v2 v1 0 |
1265  */
1266 
1267  /* assume che il Vec3 sia un'array di 3 reali */
1268  const doublereal* pdFrom = v.pGetVec();
1269 
1270  /* Coefficiente 1,2 */
1271  doublereal* pdm = pdMatm1 + iSubIt;
1272  integer* pir = piRowm1 + iSubIt;
1273  integer* pic = piColm1 + iSubIt;
1274 
1275  pdm[0] = -pdFrom[V3]; // -v.dGet(3);
1276  pir[0] = iFirstRow + 1;
1277  pic[0] = iFirstCol + 2;
1278 
1279  /* Coefficiente 1,3 */
1280  pdm[1] = pdFrom[V2]; // v.dGet(2);
1281  pir[1] = iFirstRow + 1;
1282  pic[1] = iFirstCol + 3;
1283 
1284  /* Coefficiente 2,1 */
1285  pdm[2] = pdFrom[V3]; // v.dGet(3);
1286  pir[2] = iFirstRow + 2;
1287  pic[2] = iFirstCol + 1;
1288 
1289  /* Coefficiente 2,3 */
1290  pdm[3] = -pdFrom[V1]; // -v.dGet(1);
1291  pir[3] = iFirstRow + 2;
1292  pic[3] = iFirstCol + 3;
1293 
1294  /* Coefficiente 3,1 */
1295  pdm[4] = -pdFrom[V2]; // -v.dGet(2);
1296  pir[4] = iFirstRow + 3;
1297  pic[4] = iFirstCol + 1;
1298 
1299  /* Coefficiente 3,2 */
1300  pdm[5] = pdFrom[V1]; // v.dGet(1);
1301  pir[5] = iFirstRow + 3;
1302  pic[5] = iFirstCol + 2;
1303 }
1304 
1305 
1306 void
1308 {
1309 #ifdef DEBUG
1310  IsValid();
1311 #endif /* DEBUG */
1312 
1313  ASSERT(iNumItems > 0);
1314 
1315 #if defined HAVE_MEMSET
1316  memset(pdMatm1 + 1, 0, iNumItems*sizeof(doublereal));
1317 #else /* !HAVE_MEMSET */
1318  for (integer i = iNumItems; i > 0; i--) {
1319  pdMatm1[i] = 0.;
1320  }
1321 #endif /* HAVE_MEMSET */
1322 }
1323 
1324 
1325 /* Inserisce una matrice 3x3;
1326  * si noti che non ci sono Add, Sub, ecc. perche' la filosofia
1327  * della matrice sparsa prevede che ad ogni item (riga, colonna, valore)
1328  * corrisponda un termine che poi verra' sommato ad una matrice vera,
1329  * senza controlli su eventuali duplicazioni. */
1330 void
1332  integer iFirstCol, const Mat3x3& m)
1333 {
1334 #ifdef DEBUG
1335  IsValid();
1336 
1337  ASSERT(iNumItems >= 9);
1338  ASSERT(iSubIt > 0);
1339  ASSERT(iSubIt <= iNumItems - 8);
1340  ASSERT(iFirstRow >= 0);
1341  ASSERT(iFirstCol >= 0);
1342 #endif /* DEBUG */
1343 
1344  /* Attenzione agli argomenti:
1345  * iSubIt e' il primo indice della matrice da utilizzare,
1346  * con 1 <= iSubit <= iCurSize;
1347  * iFirstRow e' il primo indice di riga -1, ovvero il
1348  * primo indice di riga della sottomatrice v/\ piena e' iFirstRow+1
1349  * iFirstCol e' il primo indice di colonna -1, ovvero il
1350  * primo indice di colonna della sottomatrice m piena e' iFirstCol+1
1351  */
1352 
1353  /* Per efficienza, vengono scritte esplicitamente tutte
1354  * le assegnazioni;
1355  * la funzione quindi non e' messa in linea intenzionalmente */
1356 
1357  /* Coefficienti 1,1-3,1 */
1358  doublereal* pdFrom = (doublereal*)m.pGetMat();
1359  doublereal* pdTmpMat = pdMatm1 + iSubIt;
1360  integer* piTmpRow = piRowm1 + iSubIt;
1361  integer* piTmpCol = piColm1 + iSubIt;
1362 
1363  iFirstRow++;
1364  iFirstCol++;
1365 
1366  /* Prima riga */
1367  pdTmpMat[0] = pdFrom[M11];
1368  piTmpRow[0] = iFirstRow++;
1369  piTmpCol[0] = iFirstCol;
1370 
1371  pdTmpMat[1] = pdFrom[M21];
1372  piTmpRow[1] = iFirstRow++;
1373  piTmpCol[1] = iFirstCol;
1374 
1375  pdTmpMat[2] = pdFrom[M31];
1376  piTmpRow[2] = iFirstRow;
1377  piTmpCol[2] = iFirstCol;
1378 
1379  /* Seconda riga */
1380  iFirstRow -= 2;
1381  iFirstCol++;
1382 
1383  pdTmpMat[3] = pdFrom[M12];
1384  piTmpRow[3] = iFirstRow++;
1385  piTmpCol[3] = iFirstCol;
1386 
1387  pdTmpMat[4] = pdFrom[M22];
1388  piTmpRow[4] = iFirstRow++;
1389  piTmpCol[4] = iFirstCol;
1390 
1391  pdTmpMat[5] = pdFrom[M32];
1392  piTmpRow[5] = iFirstRow;
1393  piTmpCol[5] = iFirstCol;
1394 
1395  /* Terza riga */
1396  iFirstRow -= 2;
1397  iFirstCol++;
1398 
1399  pdTmpMat[6] = pdFrom[M13];
1400  piTmpRow[6] = iFirstRow++;
1401  piTmpCol[6] = iFirstCol;
1402 
1403  pdTmpMat[7] = pdFrom[M23];
1404  piTmpRow[7] = iFirstRow++;
1405  piTmpCol[7] = iFirstCol;
1406 
1407  pdTmpMat[8] = pdFrom[M33];
1408  piTmpRow[8] = iFirstRow;
1409  piTmpCol[8] = iFirstCol;
1410 }
1411 
1412 
1413 /* somma la matrice ad un matrix handler usando i metodi generici */
1416 {
1417 #ifdef DEBUG
1418  IsValid();
1419  MH.IsValid();
1420 #endif /* DEBUG */
1421 
1422  for (integer i = iNumItems; i > 0; i--) {
1423  ASSERT(piRowm1[i] > 0);
1424  ASSERT(piRowm1[i] <= MH.iGetNumRows());
1425  ASSERT(piColm1[i] > 0);
1426  ASSERT(piColm1[i] <= MH.iGetNumCols());
1427 
1428  MH(piRowm1[i], piColm1[i]) += pdMatm1[i];
1429  }
1430 
1431  return MH;
1432 }
1433 
1434 
1435 /* somma la matrice, trasposta, ad un matrix handler usando i metodi generici */
1438 {
1439 #ifdef DEBUG
1440  IsValid();
1441  MH.IsValid();
1442 #endif /* DEBUG */
1443 
1444  for (integer i = iNumItems; i > 0; i--) {
1445  ASSERT(piRowm1[i] > 0);
1446  ASSERT(piRowm1[i] <= MH.iGetNumCols());
1447  ASSERT(piColm1[i] > 0);
1448  ASSERT(piColm1[i] <= MH.iGetNumRows());
1449 
1450  MH(piColm1[i], piRowm1[i]) += pdMatm1[i];
1451  }
1452 
1453  return MH;
1454 }
1455 
1456 
1457 /* somma la matrice ad un FullMatrixHandler */
1460 {
1461 #ifdef DEBUG
1462  IsValid();
1463  MH.IsValid();
1464 #endif /* DEBUG */
1465 
1466  doublereal **ppd = MH.ppdColsm1;
1467 
1468  for (integer i = iNumItems; i > 0; i--) {
1469  ASSERT(piRowm1[i] > 0);
1470  ASSERT(piRowm1[i] <= MH.iGetNumRows());
1471  ASSERT(piColm1[i] > 0);
1472  ASSERT(piColm1[i] <= MH.iGetNumCols());
1473 
1474  ppd[piColm1[i]][piRowm1[i]] += pdMatm1[i];
1475  }
1476 
1477  return MH;
1478 }
1479 
1480 
1481 /* somma la matrice, trasposta, ad un FullMatrixHandler */
1484 {
1485 #ifdef DEBUG
1486  IsValid();
1487  MH.IsValid();
1488 #endif /* DEBUG */
1489 
1490  doublereal **ppd = MH.ppdColsm1;
1491 
1492  for (integer i = iNumItems; i > 0; i--) {
1493  ASSERT(piRowm1[i] > 0);
1494  ASSERT(piRowm1[i] <= MH.iGetNumCols());
1495  ASSERT(piColm1[i] > 0);
1496  ASSERT(piColm1[i] <= MH.iGetNumRows());
1497 
1498  ppd[piRowm1[i]][piColm1[i]] += pdMatm1[i];
1499  }
1500 
1501  return MH;
1502 }
1503 
1504 
1505 /* sottrae la matrice da un matrix handler usando i metodi generici */
1508 {
1509 #ifdef DEBUG
1510  IsValid();
1511  MH.IsValid();
1512 #endif /* DEBUG */
1513 
1514  for (integer i = iNumItems; i > 0; i--) {
1515  ASSERT(piRowm1[i] > 0);
1516  ASSERT(piRowm1[i] <= MH.iGetNumRows());
1517  ASSERT(piColm1[i] > 0);
1518  ASSERT(piColm1[i] <= MH.iGetNumCols());
1519 
1520  MH(piRowm1[i], piColm1[i]) -= pdMatm1[i];
1521  }
1522 
1523  return MH;
1524 }
1525 
1526 
1527 /* sottrae la matrice, trasposta, da un matrix handler usando i metodi generici */
1530 {
1531 #ifdef DEBUG
1532  IsValid();
1533  MH.IsValid();
1534 #endif /* DEBUG */
1535 
1536  for (integer i = iNumItems; i > 0; i--) {
1537  ASSERT(piRowm1[i] > 0);
1538  ASSERT(piRowm1[i] <= MH.iGetNumCols());
1539  ASSERT(piColm1[i] > 0);
1540  ASSERT(piColm1[i] <= MH.iGetNumRows());
1541 
1542  MH(piColm1[i], piRowm1[i]) -= pdMatm1[i];
1543  }
1544 
1545  return MH;
1546 }
1547 
1548 
1549 /* sottrae la matrice da un FullMatrixHandler */
1552 {
1553 #ifdef DEBUG
1554  IsValid();
1555  MH.IsValid();
1556 #endif /* DEBUG */
1557 
1558  doublereal **ppd = MH.ppdColsm1;
1559 
1560  for (integer i = iNumItems; i > 0; i--) {
1561  ASSERT(piRowm1[i] > 0);
1562  ASSERT(piRowm1[i] <= MH.iGetNumRows());
1563  ASSERT(piColm1[i] > 0);
1564  ASSERT(piColm1[i] <= MH.iGetNumCols());
1565 
1566  ppd[piColm1[i]][piRowm1[i]] -= pdMatm1[i];
1567  }
1568 
1569  return MH;
1570 }
1571 
1572 /* sottrae la matrice, trasposta, da un FullMatrixHandler */
1575 {
1576 #ifdef DEBUG
1577  IsValid();
1578  MH.IsValid();
1579 #endif /* DEBUG */
1580 
1581  doublereal **ppd = MH.ppdColsm1;
1582 
1583  for (integer i = iNumItems; i > 0; i--) {
1584  ASSERT(piRowm1[i] > 0);
1585  ASSERT(piRowm1[i] <= MH.iGetNumCols());
1586  ASSERT(piColm1[i] > 0);
1587  ASSERT(piColm1[i] <= MH.iGetNumRows());
1588 
1589  ppd[piRowm1[i]][piColm1[i]] -= pdMatm1[i];
1590  }
1591 
1592  return MH;
1593 }
1594 
1595 /* SparseSubMatrixHandler - end */
1596 
1597 
1598 /* MySubVectorHandler - begin */
1599 
1601 : MyVectorHandler(), piRowm1(NULL) {
1602  Resize(iSize);
1603 }
1604 
1606  doublereal* pdTmpVec)
1607 : MyVectorHandler(iSize, pdTmpVec), piRowm1(0)
1608 {
1609  if (piTmpRow == 0) {
1610  /* ... because set by MyVectorHandler() */
1611  ASSERT(bOwnsMemory == true);
1612 
1613  if (pdTmpVec == 0) {
1614  silent_cerr("MySubVectorHandler(): illegal args" << std::endl);
1616  }
1617 
1618  SAFENEWARR(piRowm1, integer, iSize);
1619  }
1620 
1621  piRowm1--;
1622 #ifdef DEBUG
1623  IsValid();
1624 #endif /* DEBUG */
1625 }
1626 
1627 void
1629 {
1630  if (iSize < 0) {
1631  silent_cerr("Negative size!" << std::endl);
1633  }
1634 
1635  ASSERT((piRowm1 == NULL && pdVecm1 == NULL)
1636  || (piRowm1 != NULL && pdVecm1 != NULL));
1637 
1638  if (!bOwnsMemory && piRowm1 != NULL) {
1639  if (iSize > iMaxSize) {
1640  silent_cerr("Can't resize to " << iSize
1641  << ": larger than max size " << iMaxSize
1642  << std::endl);
1644  }
1645  iCurSize = iSize;
1646 
1647  } else {
1648  if (piRowm1 != NULL) {
1649  if (iSize <= iMaxSize) {
1650  iCurSize = iSize;
1651 
1652  } else {
1653  doublereal* pd = NULL;
1654  SAFENEWARR(pd, doublereal, iSize);
1655  pd--;
1656 
1657  integer* pi = NULL;
1658  SAFENEWARR(pi, integer, iSize);
1659  pi--;
1660 
1661  for (integer i = iCurSize; i > 0; i--) {
1662  pd[i] = pdVecm1[i];
1663  pi[i] = piRowm1[i];
1664  }
1665 
1666  pdVecm1++;
1668 
1669  piRowm1++;
1671 
1672  pdVecm1 = pd;
1673  piRowm1 = pi;
1674  iMaxSize = iCurSize = iSize;
1675  }
1676 
1677  } else {
1678  if (iSize > 0) {
1679  bOwnsMemory = true;
1680 
1681  SAFENEWARR(pdVecm1, doublereal, iSize);
1682  pdVecm1--;
1683 
1684  SAFENEWARR(piRowm1, integer, iSize);
1685  piRowm1--;
1686 
1687  iMaxSize = iCurSize = iSize;
1688  }
1689  }
1690  }
1691 }
1692 
1693 void
1695 {
1696  if (bOwnsMemory) {
1697  if (pdVecm1 != NULL) {
1698  pdVecm1++;
1700 
1701  piRowm1++;
1703  }
1704 
1705  bOwnsMemory = false;
1706  }
1707 
1708  iMaxSize = iCurSize = 0;
1709  pdVecm1 = NULL;
1710  piRowm1 = NULL;
1711 }
1712 
1713 void
1715  integer* pi, integer iMSize)
1716 {
1717  if (bOwnsMemory && pdVecm1 != NULL) {
1718  Detach();
1719  bOwnsMemory = false;
1720  }
1721 
1722  iMaxSize = iCurSize = iSize;
1723  if (iMSize >= iSize) {
1724  iMaxSize = iMSize;
1725  }
1726 
1727  pdVecm1 = pd - 1;
1728  piRowm1 = pi - 1;
1729 }
1730 
1731 #ifdef DEBUG
1732 void
1733 MySubVectorHandler::IsValid(void) const
1734 {
1735  MyVectorHandler::IsValid();
1736 
1737  ASSERT(piRowm1 != NULL);
1738 
1739 #ifdef DEBUG_MEMMANAGER
1740  ASSERT(defaultMemoryManager.fIsValid(piRowm1 + 1,
1742 #endif /* DEBUG_MEMMANAGER */
1743 }
1744 #endif /* DEBUG */
1745 
1748 {
1749 #ifdef DEBUG
1750  IsValid();
1751  VH.IsValid();
1752 #endif /* DEBUG */
1753 
1754  for (integer i = iGetSize(); i > 0; i--) {
1755 #if 0
1756  /* FIXME: workaround for SchurVectorHandler... */
1757  VH(piRowm1[i]) += pdVecm1[i];
1758 #endif
1759  VH.IncCoef(piRowm1[i], pdVecm1[i]);
1760  }
1761 
1762  return VH;
1763 }
1764 
1767 {
1768 #ifdef DEBUG
1769  IsValid();
1770  VH.IsValid();
1771 #endif /* DEBUG */
1772 
1773  doublereal* pdm1 = VH.pdGetVec() - 1;
1774  for (integer i = iGetSize(); i > 0; i--) {
1775  pdm1[piRowm1[i]] += pdVecm1[i];
1776  }
1777  return VH;
1778 }
1779 
1780 std::ostream&
1781 operator << (std::ostream& out, const SubVectorHandler& v)
1782 {
1783 #ifdef DEBUG
1784  v.IsValid();
1785 #endif /* DEBUG */
1786 
1787  integer iRow = v.iGetSize();
1788 
1789  ASSERT(iRow > 0);
1790 
1791  for (integer i = 1; i <= iRow; i++) {
1792  out << std::setw(12) << v.iGetRowIndex(i)
1793  << " " << std::setw(12) << v(i) << std::endl;
1794  }
1795 
1796  return out << std::endl;
1797 }
1798 
1799 /* MySubVectorHandler - end */
void PutT(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1338
SparseSubMatrixHandler(const SparseSubMatrixHandler &)
void PutCross(integer iFirstRow, integer iFirstCol, const Vec3 &v)
Definition: submat.cc:573
integer * piRowm1
Definition: submat.h:189
void PutMat3x3(integer iSubIt, integer iFirstRow, integer iFirstCol, const Mat3x3 &m)
Definition: submat.cc:1331
virtual integer iGetNumCols(void) const =0
std::ostream & operator<<(std::ostream &out, const FullSubMatrixHandler &m)
Definition: submat.cc:836
MatrixHandler & SubFromT(MatrixHandler &MH) const
Definition: submat.cc:748
integer * piRowm1
Definition: submat.h:768
void Attach(integer iSize, doublereal *pd, integer *pi, integer iMSize=0)
Definition: submat.cc:1714
integer iCurSize
Definition: vh.h:156
virtual ~SparseSubMatrixHandler(void)
Definition: submat.cc:1011
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
void PutT(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:239
integer iMaxCols
Definition: fullmh.h:67
void Add(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:209
void Put(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1113
Definition: matvec3.h:59
integer * piColm1
Definition: submat.h:191
void Resize(integer iNewRow, integer iNewCol)
Definition: submat.cc:138
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
virtual void IncCoef(integer iRow, const doublereal &dCoef)=0
integer iNumRows
Definition: fullmh.h:63
void Put(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:221
MatrixHandler & AddTo(MatrixHandler &MH) const
Definition: submat.cc:604
void Add(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1039
void ResizeReset(integer iNewRow, integer iNewCol)
Definition: submat.cc:1084
MatrixHandler & SubFromT(MatrixHandler &MH) const
Definition: submat.cc:1529
doublereal * pdMatm1
Definition: submat.h:772
#define NO_OP
Definition: myassert.h:74
void PutCross(integer iSubIt, integer iFirstRow, integer iFirstCol, const Vec3 &v)
Definition: submat.cc:1236
Definition: matvec3.h:50
virtual void Resize(integer iNewRows, integer iNewCols)
Definition: fullmh.cc:174
MatrixHandler & AddTo(MatrixHandler &MH) const
Definition: submat.cc:1415
virtual integer iGetSize(void) const =0
doublereal * pdVecm1
Definition: vh.h:158
void Resize(integer iNewRow, integer iNewCol)
Definition: submat.cc:1053
virtual ~FullSubMatrixHandler(void)
Definition: submat.cc:87
void AddT(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:227
Definition: matvec3.h:58
virtual void Resize(integer iSize)
Definition: submat.cc:1628
Definition: matvec3.h:51
integer iMaxSize
Definition: vh.h:155
Definition: matvec3.h:55
FullSubMatrixHandler(const FullSubMatrixHandler &)
virtual integer iGetSize(void) const
Definition: submat.h:1523
void Detach(void)
Definition: submat.cc:1694
virtual VectorHandler & AddTo(VectorHandler &VH) const
Definition: submat.cc:1747
Definition: matvec3.h:56
void Attach(int iRows, int iCols, integer *piTmpIndx)
Definition: submat.cc:193
Definition: matvec3.h:63
Definition: matvec3.h:62
Definition: matvec3.h:61
MatrixHandler & AddToT(MatrixHandler &MH) const
Definition: submat.cc:632
Definition: mbdyn.h:76
MatrixHandler & SubFrom(MatrixHandler &MH) const
Definition: submat.cc:720
integer iRawSize
Definition: fullmh.h:66
Definition: matvec3.h:57
virtual integer iGetRowIndex(integer iSubRow) const =0
void Reset(void)
Definition: submat.cc:115
integer iNumCols
Definition: fullmh.h:64
void Attach(int iNumEntr, doublereal *pdTmpMat, integer *piTmpIndx)
Definition: submat.cc:1095
virtual ~SubMatrixHandler(void)
Definition: submat.cc:43
#define ASSERT(expression)
Definition: colamd.c:977
void Sub(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1076
#define defaultMemoryManager
Definition: mynewmem.h:691
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
virtual integer iGetNumCols(void) const
Definition: fullmh.h:229
static std::stack< cleanup * > c
Definition: cleanup.cc:59
const doublereal * pGetMat(void) const
Definition: matvec3.h:743
void SubT(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:233
void PutDiag(integer iSubIt, integer iFirstRow, integer iFirstCol, const Vec3 &v)
Definition: submat.cc:1125
Definition: matvec3.h:49
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
Definition: matvec3.h:60
void PutDiag(integer iFirstRow, integer iFirstCol, const Vec3 &v)
Definition: submat.cc:523
void Sub(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:215
integer * piRowm1
Definition: submat.h:1474
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
bool bOwnsMemory
Definition: fullmh.h:61
doublereal ** ppdColsm1
Definition: fullmh.h:72
void SubT(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1301
integer * piColm1
Definition: submat.h:770
double doublereal
Definition: colamd.c:52
virtual doublereal * pdGetVec(void) const
Definition: vh.h:245
long int integer
Definition: colamd.c:51
void AddT(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1264
MatrixHandler & AddToT(MatrixHandler &MH) const
Definition: submat.cc:1437
virtual integer iGetNumRows(void) const
Definition: fullmh.h:225
virtual integer iGetNumRows(void) const =0
MySubVectorHandler(const MySubVectorHandler &)
integer iVecSize
Definition: submat.h:186
void Reset(void)
Definition: fullmh.cc:44
MatrixHandler & SubFrom(MatrixHandler &MH) const
Definition: submat.cc:1507
bool bOwnsMemory
Definition: vh.h:152