MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
fullmh.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/fullmh.cc,v 1.54 2017/01/12 14:43:53 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 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
33 
34 #include <cstring> /* for memset() */
35 #include <iostream>
36 #include <iomanip>
37 
38 #include <fullmh.h>
39 #include <submat.h>
40 
41 /* FullMatrixHandler - begin */
42 
43 void
44 FullMatrixHandler::Reset(void)
45 {
46 #ifdef DEBUG
47  IsValid();
48 #endif /* DEBUG */
49 
50  for (integer c = iNumCols; c > 0; c--) {
51  for (integer r = iNumRows; r > 0; r--) {
52  ppdColsm1[c][r] = 0.;
53  }
54  }
55 }
56 
57 void
58 FullMatrixHandler::CreateColRow(integer iNR, integer iNC)
59 {
60  ASSERT(iNC > 0);
61  ASSERT(iNC > 0);
62  ASSERT(pdRawm1);
63  ASSERT(ppdColsm1);
64 
65  doublereal *pd = pdRawm1 + iNC*iNR;
66 
67  for (integer c = iNC; c > 0; c--) {
68  pd -= iNR;
69  ppdColsm1[c] = pd;
70  }
71 
72  ASSERT(pd == pdRawm1);
73 }
74 
75 FullMatrixHandler::FullMatrixHandler(doublereal* pd, doublereal** ppd,
76  integer iSize, integer iNR, integer iNC,
77  integer iMaxC)
78 : bOwnsMemory(false),
79 iNumRows(iNR), iNumCols(iNC), iRawSize(iSize), iMaxCols(iMaxC),
80 pdRaw(pd), pdRawm1(0),
81 ppdCols(ppd), ppdColsm1(0), m_end(*this, true)
82 {
83  if (iMaxCols == 0) {
84  iMaxCols = iNumCols;
85  }
86 
87  ASSERT(pd);
88  pdRawm1 = pd - 1;
89  ASSERT(ppd);
90  ppdColsm1 = ppd - 1;
91 
92  CreateColRow(iNumRows, iNumCols);
93 
94 #ifdef DEBUG
95  IsValid();
96 #endif /* DEBUG */
97 }
98 
99 /* costruttore che si alloca la memoria */
100 FullMatrixHandler::FullMatrixHandler(integer iNR, integer iNC)
101 : bOwnsMemory(true),
102 iNumRows(iNR), iNumCols(iNC ? iNC : iNumRows),
103 iRawSize(iNumRows*iNumCols), iMaxCols(iNumCols),
104 pdRaw(NULL), pdRawm1(NULL),
105 ppdCols(NULL), ppdColsm1(NULL), m_end(*this, true)
106 {
107  ASSERT(iNumRows > 0);
108  ASSERT(iNumCols > 0);
109  Resize(iNumRows, iNumCols);
110  Reset();
111 
112 #ifdef DEBUG
113  IsValid();
114 #endif /* DEBUG */
115 }
116 
117 
118 /* costruttore che non fa nulla */
119 FullMatrixHandler::FullMatrixHandler(void)
120 : bOwnsMemory(true),
121 iNumRows(0), iNumCols(0), iRawSize(0), iMaxCols(0),
122 pdRaw(NULL), pdRawm1(NULL),
123 ppdCols(NULL), ppdColsm1(NULL), m_end(*this, true)
124 {
125  NO_OP;
126 }
127 
128 /* costruttore di copia */
129 FullMatrixHandler::FullMatrixHandler(const FullMatrixHandler& MH)
130 : bOwnsMemory(true),
131 iNumRows(0), iNumCols(0), iRawSize(0), iMaxCols(0),
132 pdRaw(NULL), pdRawm1(NULL),
133 ppdCols(NULL), ppdColsm1(NULL), m_end(*this, true)
134 {
135  if (MH.pdRaw != 0) {
136  Resize(MH.iNumRows, MH.iNumCols);
137 
138  for (integer i = 0; i < iNumRows*iNumCols; i++) {
139  pdRaw[i] = MH.pdRaw[i];
140  }
141  }
142 }
143 
144 /* Overload di = usato per l'assemblaggio delle matrici */
145 FullMatrixHandler&
146 FullMatrixHandler::operator = (const FullMatrixHandler& MH)
147 {
148  if (pdRaw == 0) {
149  Resize(MH.iNumRows, MH.iNumCols);
150 
151  } else if (iNumRows != MH.iNumRows || iNumCols != MH.iNumCols) {
152  silent_cerr("FullMatrixHandler::operator = (const FullMatrixHandler&): "
153  "incompatible size (" << iNumRows << ", " << iNumCols << ") "
154  " <> (" << MH.iNumRows << ", " << MH.iNumCols << ")"
155  << std::endl);
156  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
157  }
158 
159  for (integer i = 0; i < iNumRows*iNumCols; i++) {
160  pdRaw[i] = MH.pdRaw[i];
161  }
162 
163  return *this;
164 }
165 
166 
167 FullMatrixHandler::~FullMatrixHandler(void)
168 {
169  Detach();
170 }
171 
172 /* ridimensiona la matrice (se possiede la memoria) */
173 void
174 FullMatrixHandler::Resize(integer iNewRows, integer iNewCols)
175 {
176  if (iNewRows < 0 || iNewCols < 0) {
177  silent_cerr("Negative size!" << std::endl);
178  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
179  }
180 
181  integer iSize = iNewRows*iNewCols;
182  if (bOwnsMemory) {
183  if (pdRaw != NULL) {
184  if (iSize > iRawSize) {
185  /* crea */
186  doublereal* pd = NULL;
187  SAFENEWARR(pd, doublereal, iSize);
188  for (integer i = iRawSize; i-- > 0; ) {
189  pd[i] = pdRaw[i];
190  }
191  iRawSize = iSize;
192  SAFEDELETEARR(pdRaw);
193  pdRaw = pd;
194  pdRawm1 = pd - 1;
195  }
196 
197  if (iNewCols > iMaxCols) {
198  /* crea */
199  SAFEDELETEARR(ppdCols);
200  SAFENEWARR(ppdCols, doublereal*, iNewCols);
201  iMaxCols = iNewCols;
202  ppdColsm1 = ppdCols-1;
203  CreateColRow(iNewRows, iNewCols);
204  }
205 
206  /* aggiorna iNumCols, iNumRows */
207  if (iNewRows != iNumRows || iNewCols != iNumCols) {
208  CreateColRow(iNewRows, iNewCols);
209  }
210 
211  iNumRows = iNewRows;
212  iNumCols = iNewCols;
213 
214  } else {
215  /* crea tutto */
216  iRawSize = iSize;
217  iMaxCols = iNewCols;
218  iNumRows = iNewRows;
219  iNumCols = iNewCols;
220  SAFENEWARR(pdRaw, doublereal, iRawSize);
221  pdRawm1 = pdRaw-1;
222  SAFENEWARR(ppdCols, doublereal*, iNumCols);
223  ppdColsm1 = ppdCols-1;
224  CreateColRow(iNumRows, iNumCols);
225  }
226 
227  } else {
228  if (pdRaw != NULL) {
229  if (iSize > iRawSize || iNewCols > iMaxCols) {
230  if (iSize > iRawSize) {
231  /* errore */
232  silent_cerr("Can't resize to "
233  << iSize
234  << ": larger than max size "
235  << iRawSize << std::endl);
236  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
237 
238  } else if (iNewCols > iMaxCols) {
239  silent_cerr("Can't resize to "
240  << iNewCols
241  << " cols: larger than max "
242  "num cols "
243  << iMaxCols << std::endl);
244  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
245  }
246 
247  } else {
248  /* aggiorna */
249  if (iNewRows != iNumRows ||
250  iNewCols != iNumCols) {
251  CreateColRow(iNewRows, iNewCols);
252  }
253  iNumRows = iNewRows;
254  iNumCols = iNewCols;
255  }
256 
257  } else {
258  /* errore */
259  silent_cerr("internal error!" << std::endl);
260  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
261  }
262  }
263 }
264 
265 
266 /* si stacca dalla memoria a cui e' associato */
267 void
268 FullMatrixHandler::Detach(void)
269 {
270  if (bOwnsMemory) {
271  if (pdRaw != NULL) {
272  SAFEDELETEARR(pdRaw);
273  }
274  if (ppdCols != NULL) {
275  SAFEDELETEARR(ppdCols);
276  }
277  }
278  iRawSize = iMaxCols = iNumRows = iNumCols = 0;
279  pdRaw = pdRawm1 = NULL;
280  ppdCols = ppdColsm1 = NULL;
281 
282  /* so Resize() can be safely invoked */
283  bOwnsMemory = true;
284 }
285 
286 
287 /* Attacca un nuovo array, con n. righe, n. colonne e dim. massima;
288  * se assente, assunta = nrighe*ncolonne */
289 void
290 FullMatrixHandler::Attach(integer iNewRows, integer iNewCols,
291  doublereal* pd, doublereal** ppd,
292  integer iMSize, integer iMaxC)
293 {
294  if (bOwnsMemory || pdRaw != NULL) {
295  Detach();
296  bOwnsMemory = false;
297  }
298 
299  iNumRows = iNewRows;
300  iNumCols = iNewCols;
301  if (iMSize > 0) {
302  if (iMSize < iNumRows*iNumCols) {
303  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
304  }
305  iRawSize = iMSize;
306 
307  } else {
308  iRawSize = iNumRows*iNumCols;
309  }
310 
311  if (iMaxC > 0) {
312  if (iMaxC < iNumCols) {
313  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
314  }
315  iMaxCols = iMaxC;
316 
317  } else {
318  iMaxCols = iNumCols;
319  }
320  pdRaw = pd;
321  pdRawm1 = pdRaw - 1;
322  ppdCols = ppd;
323  ppdColsm1 = ppdCols - 1;
324 
325 #ifdef DEBUG
326  IsValid();
327 #endif /* DEBUG */
328 }
329 
330 #ifdef DEBUG
331 /* Usata per il debug */
332 void
333 FullMatrixHandler::IsValid(void) const
334 {
335  ASSERT(pdRaw != NULL);
336  ASSERT(pdRawm1 != NULL);
337  ASSERT(ppdCols != NULL);
338  ASSERT(ppdColsm1 != NULL);
339  ASSERT(iNumRows > 0);
340  ASSERT(iNumCols > 0);
341  ASSERT(iRawSize >= iNumRows*iNumCols);
342 
343 #ifdef DEBUG_MEMMANAGER
344  ASSERT(defaultMemoryManager.fIsValid(pdRaw, iRawSize*sizeof(doublereal)));
345  ASSERT(defaultMemoryManager.fIsValid(ppdCols, iNumCols*sizeof(doublereal*)));
346 #endif /* DEBUG_MEMMANAGER */
347 }
348 #endif /* DEBUG */
349 
350 
351 std::ostream&
352 operator << (std::ostream& out, const FullMatrixHandler& m)
353 {
354 #ifdef HAVE_FMTFLAGS_IN_IOS
355  std::ios::fmtflags oldbits = out.setf(std::ios::scientific);
356 #else /* !HAVE_FMTFLAGS_IN_IOS */
357  long oldbits = out.setf(ios::scientific);
358 #endif /* !HAVE_FMTFLAGS_IN_IOS */
359 
360  out << "<FullMatrixHandler> n. rows: " << m.iNumRows
361  << "; n. cols: " << m.iNumCols << std::endl;
362  for (int i = 1; i <= m.iNumRows; i++) {
363  out << "Row " << std::setw(8) << i;
364  for (int j = 1; j <= m.iNumCols; j++) {
365  out << std::setw(10) << std::setprecision(2)
366  << m.ppdColsm1[j][i];
367  }
368  out << std::endl;
369  }
370 
371  out.flags(oldbits);
372  return out;
373 }
374 
375 extern std::ostream&
376 Write(std::ostream& out,
377  const FullMatrixHandler& m,
378  const char* s,
379  const char* s2)
380 {
381 #ifdef HAVE_FMTFLAGS_IN_IOS
382  std::ios::fmtflags oldbits = out.setf(std::ios::scientific);
383 #else /* !HAVE_FMTFLAGS_IN_IOS */
384  long oldbits = out.setf(ios::scientific);
385 #endif /* !HAVE_FMTFLAGS_IN_IOS */
386 
387  if (s2 == NULL) {
388  s2 = s;
389  }
390 
391  for (int i = 1; i <= m.iNumRows; i++) {
392  for (int j = 1; j <= m.iNumCols; j++) {
393  out << std::setw(20) << std::setprecision(12)
394  << m.ppdColsm1[j][i] << s;
395  }
396  out << s2;
397  }
398 
399  out.flags(oldbits);
400  return out;
401 }
402 
403 /* Overload di += usato per l'assemblaggio delle matrici */
404 MatrixHandler&
405 FullMatrixHandler::operator += (const SubMatrixHandler& SubMH)
406 {
407  return SubMH.AddTo(*this);
408 }
409 
410 /* Overload di -= usato per l'assemblaggio delle matrici */
411 MatrixHandler&
412 FullMatrixHandler::operator -= (const SubMatrixHandler& SubMH)
413 {
414  return SubMH.SubFrom(*this);
415 }
416 
417 /* Overload di += usato per l'assemblaggio delle matrici */
418 MatrixHandler&
419 FullMatrixHandler::operator += (const VariableSubMatrixHandler& SubMH)
420 {
421  return SubMH.AddTo(*this);
422 }
423 
424 /* Overload di -= usato per l'assemblaggio delle matrici */
425 MatrixHandler&
426 FullMatrixHandler::operator -= (const VariableSubMatrixHandler& SubMH)
427 {
428  return SubMH.SubFrom(*this);
429 }
430 
431 /* Esegue il prodotto tra due matrici e se lo memorizza */
432 void
433 FullMatrixHandler::MatMul(const FullMatrixHandler& m1,
434  const FullMatrixHandler& m2)
435 {
436  if (m1.iGetNumRows() != iGetNumRows()
437  || m2.iGetNumRows() != m1.iGetNumCols()
438  || m2.iGetNumCols() != iGetNumCols())
439  {
440  silent_cerr("FullMatrixHandler::MatMul: size mismatch "
441  "this(" << iGetNumRows() << ", " << iGetNumCols() << ") "
442  "= m1(" << m1.iGetNumRows() << ", " << m1.iGetNumCols() << ") "
443  "* m2(" << m2.iGetNumRows() << ", " << m2.iGetNumCols() << ")"
444  << std::endl);
445  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
446  }
447 
448  integer nc = m1.iGetNumCols();
449 
450  for (integer iRow = 1; iRow <= iNumRows; iRow++) {
451  for (integer iCol = 1; iCol <= iNumCols; iCol++) {
452  ppdColsm1[iCol][iRow] = 0.;
453  for (integer iK = 1; iK <= nc; iK++) {
454  // ppdColsm1[iCol][iRow] += m1(iRow, iK)*m2(iK, iCol);
455  ppdColsm1[iCol][iRow] += m1.ppdColsm1[iK][iRow]*m2.ppdColsm1[iCol][iK];
456  }
457  }
458  }
459 }
460 
461 #ifdef USE_LAPACK
462 /*
463 NAME
464  DGEMM - one of the matrix-matrix operations C := alpha*op( A )*op( B ) + beta*C,
465 
466 SYNOPSIS
467  SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
468 
469  DOUBLE PRECISION ALPHA,BETA
470 
471  INTEGER K,LDA,LDB,LDC,M,N
472 
473  CHARACTER TRANSA,TRANSB
474 
475  DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
476 
477 PURPOSE
478  DGEMM performs one of the matrix-matrix operations
479 
480  where op( X ) is one of
481 
482  op( X ) = X or op( X ) = X',
483 
484  alpha and beta are scalars, and A, B and C are matrices, with op( A ) an m by k matrix, op( B ) a k by
485  n matrix and C an m by n matrix.
486 
487 ARGUMENTS
488  TRANSA - CHARACTER*1. On entry, TRANSA specifies the form of op( A ) to be used in the matrix multiplica-
489  tion as follows:
490 
491  TRANSA = 'N' or 'n', op( A ) = A.
492 
493  TRANSA = 'T' or 't', op( A ) = A'.
494 
495  TRANSA = 'C' or 'c', op( A ) = A'.
496 
497  Unchanged on exit.
498 
499  TRANSB - CHARACTER*1. On entry, TRANSB specifies the form of op( B ) to be used in the matrix multiplica-
500  tion as follows:
501 
502  TRANSB = 'N' or 'n', op( B ) = B.
503 
504  TRANSB = 'T' or 't', op( B ) = B'.
505 
506  TRANSB = 'C' or 'c', op( B ) = B'.
507 
508  Unchanged on exit.
509 
510  M - INTEGER.
511  On entry, M specifies the number of rows of the matrix op( A ) and of the matrix C. M
512  must be at least zero. Unchanged on exit.
513 
514  N - INTEGER.
515  On entry, N specifies the number of columns of the matrix op( B ) and the number of columns of
516  the matrix C. N must be at least zero. Unchanged on exit.
517 
518  K - INTEGER.
519  On entry, K specifies the number of columns of the matrix op( A ) and the number of rows of the
520  matrix op( B ). K must be at least zero. Unchanged on exit.
521 
522  ALPHA - DOUBLE PRECISION.
523  On entry, ALPHA specifies the scalar alpha. Unchanged on exit.
524 
525  A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
526  k when TRANSA = 'N' or 'n', and is m otherwise. Before entry with TRANSA = 'N' or 'n', the
527  leading m by k part of the array A must contain the matrix A, otherwise the leading k by m
528  part of the array A must contain the matrix A. Unchanged on exit.
529 
530  LDA - INTEGER.
531  On entry, LDA specifies the first dimension of A as declared in the calling (sub) program. When
532  TRANSA = 'N' or 'n' then LDA must be at least max( 1, m ), otherwise LDA must be at least max(
533  1, k ). Unchanged on exit.
534 
535  B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
536  n when TRANSB = 'N' or 'n', and is k otherwise. Before entry with TRANSB = 'N' or 'n', the
537  leading k by n part of the array B must contain the matrix B, otherwise the leading n by k
538  part of the array B must contain the matrix B. Unchanged on exit.
539 
540  LDB - INTEGER.
541  On entry, LDB specifies the first dimension of B as declared in the calling (sub) program. When
542  TRANSB = 'N' or 'n' then LDB must be at least max( 1, k ), otherwise LDB must be at least max(
543  1, n ). Unchanged on exit.
544 
545  BETA - DOUBLE PRECISION.
546  On entry, BETA specifies the scalar beta. When BETA is supplied as zero then C need not be
547  set on input. Unchanged on exit.
548 
549  C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
550  Before entry, the leading m by n part of the array C must contain the matrix C, except when
551  beta is zero, in which case C need not be set on entry. On exit, the array C is overwritten by
552  the m by n matrix ( alpha*op( A )*op( B ) + beta*C ).
553 
554  LDC - INTEGER.
555  On entry, LDC specifies the first dimension of C as declared in the calling (sub) program.
556  LDC must be at least max( 1, m ). Unchanged on exit.
557 
558  Level 3 Blas routine.
559 
560  -- Written on 8-February-1989. Jack Dongarra, Argonne National Laboratory. Iain Duff, AERE Har-
561  well. Jeremy Du Croz, Numerical Algorithms Group Ltd. Sven Hammarling, Numerical Algorithms Group
562  Ltd.
563  */
564 /*
565  Note: dgemm performs
566 
567  C := alpha*op( A )*op( B ) + beta*C,
568 
569  in MatMatMul_base and MatTMatMul_base context,
570 
571  A := *this
572  B := in
573  C := out
574 
575  op( B ) := B
576 
577  in MatMatMul_base context,
578 
579  op( A ) := A
580 
581  in MatTMatMul_base context,
582 
583  op( A ) := A^T
584 
585  out = op( A ) * B requires alpha = 1, beta = 0
586  out += op( A ) * B requires alpha = 1, beta = 1
587  out -= op( A ) * B requires alpha = -1, beta = 1
588  */
589 
590 extern "C" int
591 __FC_DECL__(dgemm)(const char *const TRANSA, const char *const TRANSB,
592  const integer *const M, const integer *const N, const integer *const K,
593  const doublereal *const ALPHA,
594  const doublereal *const A, const integer *const LDA,
595  const doublereal *const B, const integer *const LDB,
596  const doublereal *const BETA,
597  doublereal *const C, const integer *const LDC);
598 #endif // USE_LAPACK
599 
600 MatrixHandler&
601 FullMatrixHandler::MatMatMul_base(
602  void (MatrixHandler::*op)(integer iRow, integer iCol,
603  const doublereal& dCoef),
604  MatrixHandler& out, const MatrixHandler& in) const
605 {
606  const FullMatrixHandler *pin = dynamic_cast<const FullMatrixHandler *>(&in);
607  if (pin == 0) {
608  // if input matrix is not FullMatrixHandler, use generic function
609  return MatrixHandler::MatMatMul_base(op, out, in);
610  }
611 
612  integer out_nc = out.iGetNumCols();
613  integer out_nr = out.iGetNumRows();
614  integer in_nr = in.iGetNumRows();
615 
616  if (out_nr != iGetNumRows()
617  || out_nc != in.iGetNumCols()
618  || in_nr != iGetNumCols())
619  {
620  const char *strop;
621 
622  if (op == &MatrixHandler::IncCoef) {
623  strop = "+=";
624  } else if (op == &MatrixHandler::DecCoef) {
625  strop = "-=";
626  } else if (op == &MatrixHandler::PutCoef) {
627  strop = "=";
628  } else {
629  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
630  }
631 
632  silent_cerr("FullMatrixHandler::MatMatMul_base: size mismatch "
633  "out(" << out_nr << ", " << out_nc << ") "
634  << strop << " this(" << iGetNumRows() << ", " << iGetNumCols() << ") "
635  "* in(" << in_nr << ", " << in.iGetNumCols() << ")"
636  << std::endl);
637  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
638  }
639 
640  FullMatrixHandler *pout = dynamic_cast<FullMatrixHandler *>(&out);
641  if (pout == 0) {
642  // if output matrix is not FullMatrixHandler,
643  // optimize coefficient computation
644  for (integer c = 1; c <= out_nc; c++) {
645  for (integer r = 1; r <= out_nr; r++) {
646  doublereal d = 0.;
647 
648  for (integer k = 1; k <= in_nr; k++) {
649  // this(r, k) * in(k, c)
650  d += ppdColsm1[k][r]*pin->ppdColsm1[c][k];
651  }
652 
653  (out.*op)(r, c, d);
654  }
655  }
656 
657  return out;
658  }
659 
660  // if all matrices are FullMatrixHandler,
661  // either use LAPACK or optimize coefficient computation and store
662 #ifdef USE_LAPACK
663  doublereal ALPHA, BETA;
664  if (op == &MatrixHandler::IncCoef) {
665  BETA = 1.;
666  ALPHA = 1.;
667 
668  } else if (op == &MatrixHandler::DecCoef) {
669  BETA = 1.;
670  ALPHA = -1.;
671 
672  } else if (op == &MatrixHandler::PutCoef) {
673  BETA = 0.;
674  ALPHA = 1.;
675 
676  } else {
677  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
678  }
679 
680  const char *TRANSA = "N";
681  const char *TRANSB = "N";
682  __FC_DECL__(dgemm)(TRANSA, TRANSB, &out_nr, &out_nc, &in_nr,
683  &ALPHA,
684  this->pdRaw, &out_nr,
685  pin->pdRaw, &in_nr,
686  &BETA,
687  pout->pdRaw, &out_nr);
688 
689 #else // ! USE_LAPACK
690  if (op == &MatrixHandler::IncCoef) {
691  for (integer c = 1; c <= out_nc; c++) {
692  for (integer r = 1; r <= out_nr; r++) {
693  doublereal d = 0.;
694 
695  for (integer k = 1; k <= in_nr; k++) {
696  // this(r, k) * in(k, c)
697  d += ppdColsm1[k][r]*pin->ppdColsm1[c][k];
698  }
699 
700  // (out.*op)(r, c, d);
701  pout->ppdColsm1[c][r] += d;
702  }
703  }
704 
705  } else if (op == &MatrixHandler::DecCoef) {
706  for (integer c = 1; c <= out_nc; c++) {
707  for (integer r = 1; r <= out_nr; r++) {
708  doublereal d = 0.;
709 
710  for (integer k = 1; k <= in_nr; k++) {
711  // this(r, k) * in(k, c)
712  d += ppdColsm1[k][r]*pin->ppdColsm1[c][k];
713  }
714 
715  // (out.*op)(r, c, d);
716  pout->ppdColsm1[c][r] -= d;
717  }
718  }
719 
720  } else if (op == &MatrixHandler::PutCoef) {
721  for (integer c = 1; c <= out_nc; c++) {
722  for (integer r = 1; r <= out_nr; r++) {
723  doublereal d = 0.;
724 
725  for (integer k = 1; k <= in_nr; k++) {
726  // this(r, k) * in(k, c)
727  d += ppdColsm1[k][r]*pin->ppdColsm1[c][k];
728  }
729 
730  // (out.*op)(r, c, d);
731  pout->ppdColsm1[c][r] = d;
732  }
733  }
734 
735  } else {
736  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
737  }
738 #endif // ! USE_LAPACK
739 
740  return out;
741 }
742 
743 MatrixHandler&
744 FullMatrixHandler::MatTMatMul_base(
745  void (MatrixHandler::*op)(integer iRow, integer iCol,
746  const doublereal& dCoef),
747  MatrixHandler& out, const MatrixHandler& in) const
748 {
749  const FullMatrixHandler *pin = dynamic_cast<const FullMatrixHandler *>(&in);
750  if (pin == 0) {
751  // if input matrix is not FullMatrixHandler, use generic function
752  return MatrixHandler::MatTMatMul_base(op, out, in);
753  }
754 
755  integer out_nc = out.iGetNumCols();
756  integer out_nr = out.iGetNumRows();
757  integer in_nr = in.iGetNumRows();
758 
759  if (out_nr != iGetNumCols()
760  || out_nc != in.iGetNumCols()
761  || in_nr != iGetNumRows())
762  {
763  const char *strop;
764 
765  if (op == &MatrixHandler::IncCoef) {
766  strop = "+=";
767  } else if (op == &MatrixHandler::DecCoef) {
768  strop = "-=";
769  } else if (op == &MatrixHandler::PutCoef) {
770  strop = "=";
771  } else {
772  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
773  }
774 
775  silent_cerr("FullMatrixHandler::MatTMatMul_base: size mismatch "
776  "out(" << out_nr << ", " << out_nc << ") "
777  << strop << " this(" << iGetNumRows() << ", " << iGetNumCols() << ")^T "
778  "* in(" << in_nr << ", " << in.iGetNumCols() << ")"
779  << std::endl);
780  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
781  }
782 
783  FullMatrixHandler *pout = dynamic_cast<FullMatrixHandler *>(&out);
784  if (pout == 0) {
785  // if output matrix is not FullMatrixHandler,
786  // optimize coefficient computation
787  for (integer c = 1; c <= out_nc; c++) {
788  for (integer r = 1; r <= out_nr; r++) {
789  doublereal d = 0.;
790 
791  for (integer k = 1; k <= in_nr; k++) {
792  // this(k, r) * in(k, c)
793  d += ppdColsm1[r][k]*pin->ppdColsm1[c][k];
794  }
795 
796  (out.*op)(r, c, d);
797  }
798  }
799 
800  return out;
801  }
802 
803  // if all matrices are FullMatrixHandler,
804  // either use LAPACK or optimize coefficient computation and store
805 #ifdef USE_LAPACK
806  doublereal ALPHA, BETA;
807  if (op == &MatrixHandler::IncCoef) {
808  BETA = 1.;
809  ALPHA = 1.;
810 
811  } else if (op == &MatrixHandler::DecCoef) {
812  BETA = 1.;
813  ALPHA = -1.;
814 
815  } else if (op == &MatrixHandler::PutCoef) {
816  BETA = 0.;
817  ALPHA = 1.;
818 
819  } else {
820  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
821  }
822 
823  const char *TRANSA = "T";
824  const char *TRANSB = "N";
825  __FC_DECL__(dgemm)(TRANSA, TRANSB, &out_nr, &out_nc, &in_nr,
826  &ALPHA,
827  this->pdRaw, &in_nr,
828  pin->pdRaw, &in_nr,
829  &BETA,
830  pout->pdRaw, &out_nr);
831 
832 #else // ! USE_LAPACK
833  if (op == &MatrixHandler::IncCoef) {
834  for (integer c = 1; c <= out_nc; c++) {
835  for (integer r = 1; r <= out_nr; r++) {
836  doublereal d = 0.;
837 
838  for (integer k = 1; k <= in_nr; k++) {
839  // this(k, r) * in(k, c)
840  d += ppdColsm1[r][k]*pin->ppdColsm1[c][k];
841  }
842 
843  // (out.*op)(r, c, d);
844  pout->ppdColsm1[c][r] += d;
845  }
846  }
847 
848  } else if (op == &MatrixHandler::DecCoef) {
849  for (integer c = 1; c <= out_nc; c++) {
850  for (integer r = 1; r <= out_nr; r++) {
851  doublereal d = 0.;
852 
853  for (integer k = 1; k <= in_nr; k++) {
854  // this(k, r) * in(k, c)
855  d += ppdColsm1[r][k]*pin->ppdColsm1[c][k];
856  }
857 
858  // (out.*op)(r, c, d);
859  pout->ppdColsm1[c][r] -= d;
860  }
861  }
862 
863  } else if (op == &MatrixHandler::PutCoef) {
864  for (integer c = 1; c <= out_nc; c++) {
865  for (integer r = 1; r <= out_nr; r++) {
866  doublereal d = 0.;
867 
868  for (integer k = 1; k <= in_nr; k++) {
869  // this(k, r) * in(k, c)
870  d += ppdColsm1[r][k]*pin->ppdColsm1[c][k];
871  }
872 
873  // (out.*op)(r, c, d);
874  pout->ppdColsm1[c][r] = d;
875  }
876  }
877 
878  } else {
879  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
880  }
881 #endif // ! USE_LAPACK
882 
883  return out;
884 }
885 
886 VectorHandler&
887 FullMatrixHandler::MatVecMul_base(
888  void (VectorHandler::*op)(integer iRow, const doublereal& dCoef),
889  VectorHandler& out, const VectorHandler& in) const
890 {
891  integer nr = iGetNumRows();
892  integer nc = iGetNumCols();
893 
894  ASSERT(nc == in.iGetSize());
895  ASSERT(nr == out.iGetSize());
896 
897  const MyVectorHandler *pin = dynamic_cast<const MyVectorHandler *>(&in);
898  if (pin == 0) {
899  for (integer ir = 1; ir <= nr; ir++) {
900  doublereal d = 0.;
901  for (integer ic = 1; ic <= nc; ic++) {
902  // out(ir) ? this(ir, ic) * in(ic)
903  d += ppdColsm1[ic][ir]*in(ic);
904  }
905  (out.*op)(ir, d);
906  }
907 
908  return out;
909  }
910 
911  MyVectorHandler *pout = dynamic_cast<MyVectorHandler *>(&out);
912  if (pout == 0) {
913  for (integer ir = 1; ir <= nr; ir++) {
914  doublereal d = 0.;
915  for (integer ic = 1; ic <= nc; ic++) {
916  // out(ir) ? this(ir, ic) * in(ic)
917  d += ppdColsm1[ic][ir]*pin->pdVecm1[ic];
918  }
919  (out.*op)(ir, d);
920  }
921 
922  return out;
923  }
924 
925  if (op == &VectorHandler::IncCoef) {
926  for (integer ir = 1; ir <= nr; ir++) {
927  doublereal d = 0.;
928  for (integer ic = 1; ic <= nc; ic++) {
929  // out(ir) ? this(ir, ic) * in(ic)
930  d += ppdColsm1[ic][ir]*pin->pdVecm1[ic];
931  }
932  pout->pdVecm1[ir] += d;
933  }
934 
935  } else if (op == &VectorHandler::DecCoef) {
936  for (integer ir = 1; ir <= nr; ir++) {
937  doublereal d = 0.;
938  for (integer ic = 1; ic <= nc; ic++) {
939  // out(ir) ? this(ir, ic) * in(ic)
940  d += ppdColsm1[ic][ir]*pin->pdVecm1[ic];
941  }
942  pout->pdVecm1[ir] -= d;
943  }
944 
945  } else if (op == &VectorHandler::PutCoef) {
946  for (integer ir = 1; ir <= nr; ir++) {
947  doublereal d = 0.;
948  for (integer ic = 1; ic <= nc; ic++) {
949  // out(ir) ? this(ir, ic) * in(ic)
950  d += ppdColsm1[ic][ir]*pin->pdVecm1[ic];
951  }
952  pout->pdVecm1[ir] = d;
953  }
954 
955  } else {
956  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
957  }
958 
959  return out;
960 }
961 
962 VectorHandler&
963 FullMatrixHandler::MatTVecMul_base(
964  void (VectorHandler::*op)(integer iRow, const doublereal& dCoef),
965  VectorHandler& out, const VectorHandler& in) const
966 {
967  integer nr = iGetNumCols();
968  integer nc = iGetNumRows();
969 
970  ASSERT(nc == in.iGetSize());
971  ASSERT(nr == out.iGetSize());
972 
973  const MyVectorHandler *pin = dynamic_cast<const MyVectorHandler *>(&in);
974  if (pin == 0) {
975  for (integer ir = 1; ir <= nr; ir++) {
976  doublereal d = 0.;
977  for (integer ic = 1; ic <= nc; ic++) {
978  // out(ir) ? this(ic, ir) * in(ic)
979  d += ppdColsm1[ir][ic]*in(ic);
980  }
981  (out.*op)(ir, d);
982  }
983 
984  return out;
985  }
986 
987  MyVectorHandler *pout = dynamic_cast<MyVectorHandler *>(&out);
988  if (pout == 0) {
989  for (integer ir = 1; ir <= nr; ir++) {
990  doublereal d = 0.;
991  for (integer ic = 1; ic <= nc; ic++) {
992  // out(ir) ? this(ic, ir) * in(ic)
993  d += ppdColsm1[ir][ic]*pin->pdVecm1[ic];
994  }
995  (out.*op)(ir, d);
996  }
997 
998  return out;
999  }
1000 
1001  if (op == &VectorHandler::IncCoef) {
1002  for (integer ir = 1; ir <= nr; ir++) {
1003  doublereal d = 0.;
1004  for (integer ic = 1; ic <= nc; ic++) {
1005  // out(ir) ? this(ic, ir) * in(ic)
1006  d += ppdColsm1[ir][ic]*pin->pdVecm1[ic];
1007  }
1008  pout->pdVecm1[ir] += d;
1009  }
1010 
1011  } else if (op == &VectorHandler::DecCoef) {
1012  for (integer ir = 1; ir <= nr; ir++) {
1013  doublereal d = 0.;
1014  for (integer ic = 1; ic <= nc; ic++) {
1015  // out(ir) ? this(ic, ir) * in(ic)
1016  d += ppdColsm1[ir][ic]*pin->pdVecm1[ic];
1017  }
1018  pout->pdVecm1[ir] -= d;
1019  }
1020 
1021  } else if (op == &VectorHandler::PutCoef) {
1022  for (integer ir = 1; ir <= nr; ir++) {
1023  doublereal d = 0.;
1024  for (integer ic = 1; ic <= nc; ic++) {
1025  // out(ir) ? this(ic, ir) * in(ic)
1026  d += ppdColsm1[ir][ic]*pin->pdVecm1[ic];
1027  }
1028  pout->pdVecm1[ir] = d;
1029  }
1030 
1031  } else {
1032  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1033  }
1034 
1035  return out;
1036 }
1037 
1038 void
1039 FullMatrixHandler::Add(integer iRow,
1040  integer iCol,
1041  const FullMatrixHandler & source)
1042 {
1043  integer nr = source.iGetNumRows();
1044  integer nc = source.iGetNumCols();
1045  iRow--;
1046  iCol--;
1047 
1048 #ifdef DEBUG
1049  IsValid();
1050 
1051  ASSERT(iRow >= 0);
1052  ASSERT(iCol >= 0);
1053  ASSERT(iRow + nr <= iGetNumRows());
1054  ASSERT(iCol + nc <= iGetNumCols());
1055 #endif // DEBUG
1056 #if 0
1057  if (iRow < 0 || iCol < 0
1058  || iRow + nr > iGetNumRows()
1059  || iCol + nc > iGetNumCols())
1060  {
1061  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1062  }
1063 #endif
1064 
1065  for (integer ir = 1; ir <= nr; ir++) {
1066  for (integer ic = 1; ic <= nc; ic++) {
1067  // IncCoef(ir + iRow, ic + iCol, source(ir, ic) * dCoef);
1068  ppdColsm1[ic + iCol][ir + iRow] += source.ppdColsm1[ic][ir];
1069  }
1070  }
1071 
1072  return;
1073 }
1074 
1075 void
1076 FullMatrixHandler::Sub(integer iRow,
1077  integer iCol,
1078  const FullMatrixHandler & source)
1079 {
1080  integer nr = source.iGetNumRows();
1081  integer nc = source.iGetNumCols();
1082  iRow--;
1083  iCol--;
1084 
1085 #ifdef DEBUG
1086  IsValid();
1087 
1088  ASSERT(iRow >= 0);
1089  ASSERT(iCol >= 0);
1090  ASSERT(iRow + nr <= iGetNumRows());
1091  ASSERT(iCol + nc <= iGetNumCols());
1092 #endif // DEBUG
1093 #if 0
1094  if (iRow < 0 || iCol < 0
1095  || iRow + nr > iGetNumRows()
1096  || iCol + nc > iGetNumCols())
1097  {
1098  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1099  }
1100 #endif
1101 
1102  for (integer ir = 1; ir <= nr; ir++) {
1103  for (integer ic = 1; ic <= nc; ic++) {
1104  // IncCoef(ir + iRow, ic + iCol, source(ir, ic) * dCoef);
1105  ppdColsm1[ic + iCol][ir + iRow] -= source.ppdColsm1[ic][ir];
1106  }
1107  }
1108 
1109  return;
1110 }
1111 
1112 void
1113 FullMatrixHandler::Put(integer iRow,
1114  integer iCol,
1115  const FullMatrixHandler & source)
1116 {
1117  integer nr = source.iGetNumRows();
1118  integer nc = source.iGetNumCols();
1119  iRow--;
1120  iCol--;
1121 
1122 #ifdef DEBUG
1123  IsValid();
1124 
1125  ASSERT(iRow >= 0);
1126  ASSERT(iCol >= 0);
1127  ASSERT(iRow + nr <= iGetNumRows());
1128  ASSERT(iCol + nc <= iGetNumCols());
1129 #endif // DEBUG
1130 #if 0
1131  if (iRow < 0 || iCol < 0
1132  || iRow + nr > iGetNumRows()
1133  || iCol + nc > iGetNumCols())
1134  {
1135  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1136  }
1137 #endif
1138 
1139  for (integer ir = 1; ir <= nr; ir++) {
1140  for (integer ic = 1; ic <= nc; ic++) {
1141  // IncCoef(ir + iRow, ic + iCol, source(ir, ic) * dCoef);
1142  ppdColsm1[ic + iCol][ir + iRow] = source.ppdColsm1[ic][ir];
1143  }
1144  }
1145 
1146  return;
1147 }
1148 
1149 void
1150 FullMatrixHandler::Add(integer iRow,
1151  integer iCol,
1152  const FullMatrixHandler & source,
1153  const doublereal dCoef)
1154 {
1155  integer nr = source.iGetNumRows();
1156  integer nc = source.iGetNumCols();
1157  iRow--;
1158  iCol--;
1159 
1160 #ifdef DEBUG
1161  IsValid();
1162 
1163  ASSERT(iRow >= 0);
1164  ASSERT(iCol >= 0);
1165  ASSERT(iRow + nr <= iGetNumRows());
1166  ASSERT(iCol + nc <= iGetNumCols());
1167 #endif // DEBUG
1168 #if 0
1169  if (iRow < 0 || iCol < 0
1170  || iRow + nr > iGetNumRows()
1171  || iCol + nc > iGetNumCols())
1172  {
1173  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1174  }
1175 #endif
1176 
1177  for (integer ir = 1; ir <= nr; ir++) {
1178  for (integer ic = 1; ic <= nc; ic++) {
1179  // IncCoef(ir + iRow, ic + iCol, source(ir, ic) * dCoef);
1180  ppdColsm1[ic + iCol][ir + iRow] += dCoef*source.ppdColsm1[ic][ir];
1181  }
1182  }
1183 
1184  return;
1185 }
1186 
1187 void
1188 FullMatrixHandler::Sub(integer iRow,
1189  integer iCol,
1190  const FullMatrixHandler & source,
1191  const doublereal dCoef)
1192 {
1193  integer nr = source.iGetNumRows();
1194  integer nc = source.iGetNumCols();
1195  iRow--;
1196  iCol--;
1197 
1198 #ifdef DEBUG
1199  IsValid();
1200 
1201  ASSERT(iRow >= 0);
1202  ASSERT(iCol >= 0);
1203  ASSERT(iRow + nr <= iGetNumRows());
1204  ASSERT(iCol + nc <= iGetNumCols());
1205 #endif // DEBUG
1206 #if 0
1207  if (iRow < 0 || iCol < 0
1208  || iRow + nr > iGetNumRows()
1209  || iCol + nc > iGetNumCols())
1210  {
1211  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1212  }
1213 #endif
1214 
1215  for (integer ir = 1; ir <= nr; ir++) {
1216  for (integer ic = 1; ic <= nc; ic++) {
1217  // IncCoef(ir + iRow, ic + iCol, source(ir, ic) * dCoef);
1218  ppdColsm1[ic + iCol][ir + iRow] -= dCoef*source.ppdColsm1[ic][ir];
1219  }
1220  }
1221 
1222  return;
1223 }
1224 
1225 void
1226 FullMatrixHandler::Put(integer iRow,
1227  integer iCol,
1228  const FullMatrixHandler & source,
1229  const doublereal dCoef)
1230 {
1231  integer nr = source.iGetNumRows();
1232  integer nc = source.iGetNumCols();
1233  iRow--;
1234  iCol--;
1235 
1236 #ifdef DEBUG
1237  IsValid();
1238 
1239  ASSERT(iRow >= 0);
1240  ASSERT(iCol >= 0);
1241  ASSERT(iRow + nr <= iGetNumRows());
1242  ASSERT(iCol + nc <= iGetNumCols());
1243 #endif // DEBUG
1244 #if 0
1245  if (iRow < 0 || iCol < 0
1246  || iRow + nr > iGetNumRows()
1247  || iCol + nc > iGetNumCols())
1248  {
1249  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1250  }
1251 #endif
1252 
1253  for (integer ir = 1; ir <= nr; ir++) {
1254  for (integer ic = 1; ic <= nc; ic++) {
1255  // IncCoef(ir + iRow, ic + iCol, source(ir, ic) * dCoef);
1256  ppdColsm1[ic + iCol][ir + iRow] = dCoef*source.ppdColsm1[ic][ir];
1257  }
1258  }
1259 
1260  return;
1261 }
1262 
1263 void
1264 FullMatrixHandler::AddT(integer iRow,
1265  integer iCol,
1266  const FullMatrixHandler & source)
1267 {
1268  integer nr = source.iGetNumCols();
1269  integer nc = source.iGetNumRows();
1270  iRow--;
1271  iCol--;
1272 
1273 #ifdef DEBUG
1274  IsValid();
1275 
1276  ASSERT(iRow >= 0);
1277  ASSERT(iCol >= 0);
1278  ASSERT(iRow + nr <= iGetNumRows());
1279  ASSERT(iCol + nc <= iGetNumCols());
1280 #endif // DEBUG
1281 #if 0
1282  if (iRow < 0 || iCol < 0
1283  || iRow + nr > iGetNumRows()
1284  || iCol + nc > iGetNumCols())
1285  {
1286  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1287  }
1288 #endif
1289 
1290  for (integer ir = 1; ir <= nr; ir++) {
1291  for (integer ic = 1; ic <= nc; ic++) {
1292  // IncCoef(ir + iRow, ic + iCol, source(ic, ir) * dCoef);
1293  ppdColsm1[ic + iCol][ir + iRow] += source.ppdColsm1[ir][ic];
1294  }
1295  }
1296 
1297  return;
1298 }
1299 
1300 void
1301 FullMatrixHandler::SubT(integer iRow,
1302  integer iCol,
1303  const FullMatrixHandler & source)
1304 {
1305  integer nr = source.iGetNumCols();
1306  integer nc = source.iGetNumRows();
1307  iRow--;
1308  iCol--;
1309 
1310 #ifdef DEBUG
1311  IsValid();
1312 
1313  ASSERT(iRow >= 0);
1314  ASSERT(iCol >= 0);
1315  ASSERT(iRow + nr <= iGetNumRows());
1316  ASSERT(iCol + nc <= iGetNumCols());
1317 #endif // DEBUG
1318 #if 0
1319  if (iRow < 0 || iCol < 0
1320  || iRow + nr > iGetNumRows()
1321  || iCol + nc > iGetNumCols())
1322  {
1323  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1324  }
1325 #endif
1326 
1327  for (integer ir = 1; ir <= nr; ir++) {
1328  for (integer ic = 1; ic <= nc; ic++) {
1329  // IncCoef(ir + iRow, ic + iCol, source(ic, ir) * dCoef);
1330  ppdColsm1[ic + iCol][ir + iRow] -= source.ppdColsm1[ir][ic];
1331  }
1332  }
1333 
1334  return;
1335 }
1336 
1337 void
1338 FullMatrixHandler::PutT(integer iRow,
1339  integer iCol,
1340  const FullMatrixHandler & source)
1341 {
1342  integer nr = source.iGetNumCols();
1343  integer nc = source.iGetNumRows();
1344  iRow--;
1345  iCol--;
1346 
1347 #ifdef DEBUG
1348  IsValid();
1349 
1350  ASSERT(iRow >= 0);
1351  ASSERT(iCol >= 0);
1352  ASSERT(iRow + nr <= iGetNumRows());
1353  ASSERT(iCol + nc <= iGetNumCols());
1354 #endif // DEBUG
1355 #if 0
1356  if (iRow < 0 || iCol < 0
1357  || iRow + nr > iGetNumRows()
1358  || iCol + nc > iGetNumCols())
1359  {
1360  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1361  }
1362 #endif
1363 
1364  for (integer ir = 1; ir <= nr; ir++) {
1365  for (integer ic = 1; ic <= nc; ic++) {
1366  // IncCoef(ir + iRow, ic + iCol, source(ic, ir) * dCoef);
1367  ppdColsm1[ic + iCol][ir + iRow] = source.ppdColsm1[ir][ic];
1368  }
1369  }
1370 
1371  return;
1372 }
1373 
1374 void
1375 FullMatrixHandler::AddT(integer iRow,
1376  integer iCol,
1377  const FullMatrixHandler & source,
1378  const doublereal dCoef)
1379 {
1380  integer nr = source.iGetNumCols();
1381  integer nc = source.iGetNumRows();
1382  iRow--;
1383  iCol--;
1384 
1385 #ifdef DEBUG
1386  IsValid();
1387 
1388  ASSERT(iRow >= 0);
1389  ASSERT(iCol >= 0);
1390  ASSERT(iRow + nr <= iGetNumRows());
1391  ASSERT(iCol + nc <= iGetNumCols());
1392 #endif // DEBUG
1393 #if 0
1394  if (iRow < 0 || iCol < 0
1395  || iRow + nr > iGetNumRows()
1396  || iCol + nc > iGetNumCols())
1397  {
1398  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1399  }
1400 #endif
1401 
1402  for (integer ir = 1; ir <= nr; ir++) {
1403  for (integer ic = 1; ic <= nc; ic++) {
1404  // IncCoef(ir + iRow, ic + iCol, source(ic, ir) * dCoef);
1405  ppdColsm1[ic + iCol][ir + iRow] += dCoef*source.ppdColsm1[ir][ic];
1406  }
1407  }
1408 
1409  return;
1410 }
1411 
1412 void
1413 FullMatrixHandler::SubT(integer iRow,
1414  integer iCol,
1415  const FullMatrixHandler & source,
1416  const doublereal dCoef)
1417 {
1418  integer nr = source.iGetNumCols();
1419  integer nc = source.iGetNumRows();
1420  iRow--;
1421  iCol--;
1422 
1423 #ifdef DEBUG
1424  IsValid();
1425 
1426  ASSERT(iRow >= 0);
1427  ASSERT(iCol >= 0);
1428  ASSERT(iRow + nr <= iGetNumRows());
1429  ASSERT(iCol + nc <= iGetNumCols());
1430 #endif // DEBUG
1431 #if 0
1432  if (iRow < 0 || iCol < 0
1433  || iRow + nr > iGetNumRows()
1434  || iCol + nc > iGetNumCols())
1435  {
1436  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1437  }
1438 #endif
1439 
1440  for (integer ir = 1; ir <= nr; ir++) {
1441  for (integer ic = 1; ic <= nc; ic++) {
1442  // IncCoef(ir + iRow, ic + iCol, source(ic, ir) * dCoef);
1443  ppdColsm1[ic + iCol][ir + iRow] -= dCoef*source.ppdColsm1[ir][ic];
1444  }
1445  }
1446 
1447  return;
1448 }
1449 
1450 void
1451 FullMatrixHandler::PutT(integer iRow,
1452  integer iCol,
1453  const FullMatrixHandler & source,
1454  const doublereal dCoef)
1455 {
1456  integer nr = source.iGetNumCols();
1457  integer nc = source.iGetNumRows();
1458  iRow--;
1459  iCol--;
1460 
1461 #ifdef DEBUG
1462  IsValid();
1463 
1464  ASSERT(iRow >= 0);
1465  ASSERT(iCol >= 0);
1466  ASSERT(iRow + nr <= iGetNumRows());
1467  ASSERT(iCol + nc <= iGetNumCols());
1468 #endif // DEBUG
1469 #if 0
1470  if (iRow < 0 || iCol < 0
1471  || iRow + nr > iGetNumRows()
1472  || iCol + nc > iGetNumCols())
1473  {
1474  throw ErrGeneric(MBDYN_EXCEPT_ARGS);
1475  }
1476 #endif
1477 
1478  for (integer ir = 1; ir <= nr; ir++) {
1479  for (integer ic = 1; ic <= nc; ic++) {
1480  // IncCoef(ir + iRow, ic + iCol, source(ic, ir) * dCoef);
1481  ppdColsm1[ic + iCol][ir + iRow] = dCoef*source.ppdColsm1[ir][ic];
1482  }
1483  }
1484 
1485  return;
1486 }
1487 
1488 /* somma un vettore di tipo Vec3 in una data posizione */
1489 void
1490 FullMatrixHandler::Add(integer iRow, integer iCol, const Vec3& v)
1491 {
1492  /* iRow e iCol sono gli indici effettivi di riga e colonna
1493  * es. per il primo coefficiente:
1494  * iRow = 1, iCol = 1 */
1495 
1496 #ifdef DEBUG
1497  IsValid();
1498 
1499  ASSERT(iRow > 0);
1500  ASSERT(iRow <= iNumRows - 2);
1501  ASSERT(iCol > 0);
1502  ASSERT(iCol <= iNumCols);
1503 #endif /* DEBUG */
1504 
1505  /* Nota: assume che la matrice sia organizzata
1506  * "per colonne" (stile FORTRAN)
1507  */
1508  const doublereal* pdFrom = v.pGetVec();
1509 
1510  ppdColsm1[iCol][iRow] += pdFrom[V1];
1511 
1512  ppdColsm1[iCol][iRow + 1] += pdFrom[V2];
1513 
1514  ppdColsm1[iCol][iRow + 2] += pdFrom[V3];
1515 }
1516 
1517 /* sottrae un vettore di tipo Vec3 in una data posizione */
1518 void
1519 FullMatrixHandler::Sub(integer iRow, integer iCol, const Vec3& v)
1520 {
1521  /* iRow e iCol sono gli indici effettivi di riga e colonna
1522  * es. per il primo coefficiente:
1523  * iRow = 1, iCol = 1 */
1524 
1525 #ifdef DEBUG
1526  IsValid();
1527 
1528  ASSERT(iRow > 0);
1529  ASSERT(iRow <= iNumRows - 2);
1530  ASSERT(iCol > 0);
1531  ASSERT(iCol <= iNumCols);
1532 #endif /* DEBUG */
1533 
1534  /* Nota: assume che la matrice sia organizzata
1535  * "per colonne" (stile FORTRAN)
1536  */
1537  const doublereal* pdFrom = v.pGetVec();
1538 
1539  ppdColsm1[iCol][iRow] -= pdFrom[V1];
1540 
1541  ppdColsm1[iCol][iRow + 1] -= pdFrom[V2];
1542 
1543  ppdColsm1[iCol][iRow + 2] -= pdFrom[V3];
1544 }
1545 
1546 /* scrive un vettore di tipo Vec3 in una data posizione */
1547 void
1548 FullMatrixHandler::Put(integer iRow, integer iCol, const Vec3& v)
1549 {
1550  /* iRow e iCol sono gli indici effettivi di riga e colonna
1551  * es. per il primo coefficiente:
1552  * iRow = 1, iCol = 1 */
1553 
1554 #ifdef DEBUG
1555  IsValid();
1556 
1557  ASSERT(iRow > 0);
1558  ASSERT(iRow <= iNumRows - 2);
1559  ASSERT(iCol > 0);
1560  ASSERT(iCol <= iNumCols);
1561 #endif /* DEBUG */
1562 
1563  /* Nota: assume che la matrice sia organizzata
1564  * "per colonne" (stile FORTRAN)
1565  */
1566  const doublereal* pdFrom = v.pGetVec();
1567 
1568  ppdColsm1[iCol][iRow] = pdFrom[V1];
1569 
1570  ppdColsm1[iCol][iRow + 1] = pdFrom[V2];
1571 
1572  ppdColsm1[iCol][iRow + 2] = pdFrom[V3];
1573 }
1574 
1575 /* somma un vettore di tipo Vec3 trasposto in una data posizione */
1576 void
1577 FullMatrixHandler::AddT(integer iRow, integer iCol, const Vec3& v)
1578 {
1579  /* iRow e iCol sono gli indici effettivi di riga e colonna
1580  * es. per il primo coefficiente:
1581  * iRow = 1, iCol = 1 */
1582 
1583 #ifdef DEBUG
1584  IsValid();
1585 
1586  ASSERT(iRow > 0);
1587  ASSERT(iRow <= iNumRows);
1588  ASSERT(iCol > 0);
1589  ASSERT(iCol <= iNumCols - 2);
1590 #endif /* DEBUG */
1591 
1592  /* Nota: assume che la matrice sia organizzata
1593  * "per colonne" (stile FORTRAN)
1594  */
1595  const doublereal* pdFrom = v.pGetVec();
1596 
1597  ppdColsm1[iCol][iRow] += pdFrom[V1];
1598 
1599  ppdColsm1[iCol + 1][iRow] += pdFrom[V2];
1600 
1601  ppdColsm1[iCol + 2][iRow] += pdFrom[V3];
1602 }
1603 
1604 /* sottrae un vettore di tipo Vec3 trasposto in una data posizione */
1605 void
1606 FullMatrixHandler::SubT(integer iRow, integer iCol, const Vec3& v)
1607 {
1608  /* iRow e iCol sono gli indici effettivi di riga e colonna
1609  * es. per il primo coefficiente:
1610  * iRow = 1, iCol = 1 */
1611 
1612 #ifdef DEBUG
1613  IsValid();
1614 
1615  ASSERT(iRow > 0);
1616  ASSERT(iRow <= iNumRows);
1617  ASSERT(iCol > 0);
1618  ASSERT(iCol <= iNumCols - 2);
1619 #endif /* DEBUG */
1620 
1621  /* Nota: assume che la matrice sia organizzata
1622  * "per colonne" (stile FORTRAN)
1623  */
1624  const doublereal* pdFrom = v.pGetVec();
1625 
1626  ppdColsm1[iCol][iRow] -= pdFrom[V1];
1627 
1628  ppdColsm1[iCol + 1][iRow] -= pdFrom[V2];
1629 
1630  ppdColsm1[iCol + 2][iRow] -= pdFrom[V3];
1631 }
1632 
1633 /* scrive un vettore di tipo Vec3 trasposto in una data posizione */
1634 void
1635 FullMatrixHandler::PutT(integer iRow, integer iCol, const Vec3& v)
1636 {
1637  /* iRow e iCol sono gli indici effettivi di riga e colonna
1638  * es. per il primo coefficiente:
1639  * iRow = 1, iCol = 1 */
1640 
1641 #ifdef DEBUG
1642  IsValid();
1643 
1644  ASSERT(iRow > 0);
1645  ASSERT(iRow <= iNumRows);
1646  ASSERT(iCol > 0);
1647  ASSERT(iCol <= iNumCols - 2);
1648 #endif /* DEBUG */
1649 
1650  /* Nota: assume che la matrice sia organizzata
1651  * "per colonne" (stile FORTRAN)
1652  */
1653  const doublereal* pdFrom = v.pGetVec();
1654 
1655  ppdColsm1[iCol][iRow] = pdFrom[V1];
1656 
1657  ppdColsm1[iCol + 1][iRow] = pdFrom[V2];
1658 
1659  ppdColsm1[iCol + 2][iRow] = pdFrom[V3];
1660 }
1661 
1662 #if 0 /* FIXME: replace original? */
1663 /* somma un vettore di tipo Vec3 in una data posizione in diagonale */
1664 void
1665 FullMatrixHandler::AddDiag(integer iRow, integer iCol, const Vec3& v)
1666 {
1667  /* iRow e iCol sono gli indici effettivi di riga e colonna
1668  * es. per il primo coefficiente:
1669  * iRow = 1, iCol = 1 */
1670 
1671 #ifdef DEBUG
1672  IsValid();
1673 
1674  ASSERT(iRow > 0);
1675  ASSERT(iRow <= iNumRows - 2);
1676  ASSERT(iCol > 0);
1677  ASSERT(iCol <= iNumCols - 2);
1678 #endif /* DEBUG */
1679 
1680  /* Nota: assume che la matrice sia organizzata
1681  * "per colonne" (stile FORTRAN)
1682  */
1683  const doublereal* pdFrom = v.pGetVec();
1684 
1685  ppdColsm1[iCol][iRow] += pdFrom[V1];
1686 
1687  ppdColsm1[iCol + 1][iRow + 1] += pdFrom[V2];
1688 
1689  ppdColsm1[iCol + 2][iRow + 2] += pdFrom[V3];
1690 }
1691 
1692 /* sottrae un vettore di tipo Vec3 in una data posizione in diagonale */
1693 void
1694 FullMatrixHandler::SubDiag(integer iRow, integer iCol, const Vec3& v)
1695 {
1696  /* iRow e iCol sono gli indici effettivi di riga e colonna
1697  * es. per il primo coefficiente:
1698  * iRow = 1, iCol = 1 */
1699 
1700 #ifdef DEBUG
1701  IsValid();
1702 
1703  ASSERT(iRow > 0);
1704  ASSERT(iRow <= iNumRows - 2);
1705  ASSERT(iCol > 0);
1706  ASSERT(iCol <= iNumCols - 2);
1707 #endif /* DEBUG */
1708 
1709  /* Nota: assume che la matrice sia organizzata
1710  * "per colonne" (stile FORTRAN)
1711  */
1712  const doublereal* pdFrom = v.pGetVec();
1713 
1714  ppdColsm1[iCol][iRow] -= pdFrom[V1];
1715 
1716  ppdColsm1[iCol + 1][iRow + 1] -= pdFrom[V2];
1717 
1718  ppdColsm1[iCol + 2][iRow + 2] -= pdFrom[V3];
1719 }
1720 
1721 /* scrive un vettore di tipo Vec3 in una data posizione in diagonale */
1722 void
1723 FullMatrixHandler::PutDiag(integer iRow, integer iCol, const Vec3& v)
1724 {
1725  /* iRow e iCol sono gli indici effettivi di riga e colonna
1726  * es. per il primo coefficiente:
1727  * iRow = 1, iCol = 1 */
1728 
1729 #ifdef DEBUG
1730  IsValid();
1731 
1732  ASSERT(iRow > 0);
1733  ASSERT(iRow <= iNumRows - 2);
1734  ASSERT(iCol > 0);
1735  ASSERT(iCol <= iNumCols - 2);
1736 #endif /* DEBUG */
1737 
1738  /* Nota: assume che la matrice sia organizzata
1739  * "per colonne" (stile FORTRAN)
1740  */
1741  const doublereal* pdFrom = v.pGetVec();
1742 
1743  ppdColsm1[iCol][iRow] = pdFrom[V1];
1744 
1745  ppdColsm1[iCol + 1][iRow + 1] = pdFrom[V2];
1746 
1747  ppdColsm1[iCol + 2][iRow + 2] = pdFrom[V3];
1748 }
1749 
1750 /* somma un vettore di tipo Vec3 in una data posizione [ v x ] */
1751 void
1752 FullMatrixHandler::AddCross(integer iRow, integer iCol, const Vec3& v)
1753 {
1754  /* iRow e iCol sono gli indici effettivi di riga e colonna
1755  * es. per il primo coefficiente:
1756  * iRow = 1, iCol = 1 */
1757 
1758 #ifdef DEBUG
1759  IsValid();
1760 
1761  ASSERT(iRow > 0);
1762  ASSERT(iRow <= iNumRows - 2);
1763  ASSERT(iCol > 0);
1764  ASSERT(iCol <= iNumCols - 2);
1765 #endif /* DEBUG */
1766 
1767  /* Nota: assume che la matrice sia organizzata
1768  * "per colonne" (stile FORTRAN)
1769  */
1770  const doublereal* pdFrom = v.pGetVec();
1771 
1772  ppdColsm1[iCol + 1][iRow] -= pdFrom[V3];
1773  ppdColsm1[iCol + 2][iRow] += pdFrom[V2];
1774 
1775  ppdColsm1[iCol][iRow + 1] += pdFrom[V3];
1776  ppdColsm1[iCol + 2][iRow + 1] -= pdFrom[V1];
1777 
1778  ppdColsm1[iCol][iRow + 2] -= pdFrom[V2];
1779  ppdColsm1[iCol + 1][iRow + 2] += pdFrom[V1];
1780 }
1781 
1782 /* sottrae un vettore di tipo Vec3 in una data posizione [ v x ] */
1783 void
1784 FullMatrixHandler::SubCross(integer iRow, integer iCol, const Vec3& v)
1785 {
1786  /* iRow e iCol sono gli indici effettivi di riga e colonna
1787  * es. per il primo coefficiente:
1788  * iRow = 1, iCol = 1 */
1789 
1790 #ifdef DEBUG
1791  IsValid();
1792 
1793  ASSERT(iRow > 0);
1794  ASSERT(iRow <= iNumRows - 2);
1795  ASSERT(iCol > 0);
1796  ASSERT(iCol <= iNumCols - 2);
1797 #endif /* DEBUG */
1798 
1799  /* Nota: assume che la matrice sia organizzata
1800  * "per colonne" (stile FORTRAN)
1801  */
1802  const doublereal* pdFrom = v.pGetVec();
1803 
1804  ppdColsm1[iCol + 1][iRow] += pdFrom[V3];
1805  ppdColsm1[iCol + 2][iRow] -= pdFrom[V2];
1806 
1807  ppdColsm1[iCol][iRow + 1] -= pdFrom[V3];
1808  ppdColsm1[iCol + 2][iRow + 1] += pdFrom[V1];
1809 
1810  ppdColsm1[iCol][iRow + 2] += pdFrom[V2];
1811  ppdColsm1[iCol + 1][iRow + 2] -= pdFrom[V1];
1812 }
1813 
1814 /* scrive un vettore di tipo Vec3 in una data posizione [ v x ] */
1815 void
1816 FullMatrixHandler::PutCross(integer iRow, integer iCol, const Vec3& v)
1817 {
1818  /* iRow e iCol sono gli indici effettivi di riga e colonna
1819  * es. per il primo coefficiente:
1820  * iRow = 1, iCol = 1 */
1821 
1822 #ifdef DEBUG
1823  IsValid();
1824 
1825  ASSERT(iRow > 0);
1826  ASSERT(iRow <= iNumRows - 2);
1827  ASSERT(iCol > 0);
1828  ASSERT(iCol <= iNumCols - 2);
1829 #endif /* DEBUG */
1830 
1831  /* Nota: assume che la matrice sia organizzata
1832  * "per colonne" (stile FORTRAN)
1833  */
1834  const doublereal* pdFrom = v.pGetVec();
1835 
1836  ppdColsm1[iCol + 1][iRow] = -pdFrom[V3];
1837  ppdColsm1[iCol + 2][iRow] = pdFrom[V2];
1838 
1839  ppdColsm1[iCol][iRow + 1] = pdFrom[V3];
1840  ppdColsm1[iCol + 2][iRow + 1] = -pdFrom[V1];
1841 
1842  ppdColsm1[iCol][iRow + 2] = -pdFrom[V2];
1843  ppdColsm1[iCol + 1][iRow + 2] = pdFrom[V1];
1844 }
1845 #endif
1846 
1847 /* somma una matrice di tipo Mat3x3 in una data posizione */
1848 void
1849 FullMatrixHandler::Add(integer iRow, integer iCol, const Mat3x3& m)
1850 {
1851  /* iRow e iCol sono gli indici effettivi di riga e colonna
1852  * es. per il primo coefficiente:
1853  * iRow = 1, iCol = 1 */
1854 
1855 #ifdef DEBUG
1856  IsValid();
1857 
1858  ASSERT(iRow > 0);
1859  ASSERT(iRow <= iNumRows - 2);
1860  ASSERT(iCol > 0);
1861  ASSERT(iCol <= iNumCols - 2);
1862 #endif /* DEBUG */
1863 
1864  /* Nota: assume che la matrice sia organizzata
1865  * "per colonne" (stile FORTRAN) e assume
1866  * che la matrice Mat3x3 sia organizzata anch'essa per colonne */
1867  const doublereal* pdFrom = m.pGetMat();
1868 
1869  ppdColsm1[iCol][iRow] += pdFrom[M11];
1870  ppdColsm1[iCol][iRow + 1] += pdFrom[M21];
1871  ppdColsm1[iCol][iRow + 2] += pdFrom[M31];
1872 
1873  ppdColsm1[iCol + 1][iRow] += pdFrom[M12];
1874  ppdColsm1[iCol + 1][iRow + 1] += pdFrom[M22];
1875  ppdColsm1[iCol + 1][iRow + 2] += pdFrom[M32];
1876 
1877  ppdColsm1[iCol + 2][iRow] += pdFrom[M13];
1878  ppdColsm1[iCol + 2][iRow + 1] += pdFrom[M23];
1879  ppdColsm1[iCol + 2][iRow + 2] += pdFrom[M33];
1880 
1881 #if 0
1882  for (unsigned c = 0; c < 2; c++) {
1883  for (unsigned r = 0; r < 2; r++) {
1884  ppdColsm1[iCol + c][iRow + r] += pdFrom[r];
1885  }
1886  pdFrom += 3;
1887  }
1888 #endif
1889 }
1890 
1891 void
1892 FullMatrixHandler::AddT(integer iRow, integer iCol, const Mat3x3& m)
1893 {
1894  /* iRow e iCol sono gli indici effettivi di riga e colonna
1895  * es. per il primo coefficiente:
1896  * iRow = 1, iCol = 1 */
1897 
1898 #ifdef DEBUG
1899  IsValid();
1900 
1901  ASSERT(iRow > 0);
1902  ASSERT(iRow <= iNumRows - 2);
1903  ASSERT(iCol > 0);
1904  ASSERT(iCol <= iNumCols - 2);
1905 #endif /* DEBUG */
1906 
1907  /* Nota: assume che la matrice sia organizzata
1908  * "per colonne" (stile FORTRAN) e assume
1909  * che la matrice Mat3x3 sia organizzata anch'essa per colonne */
1910  const doublereal* pdFrom = m.pGetMat();
1911 
1912  ppdColsm1[iCol][iRow] += pdFrom[M11];
1913  ppdColsm1[iCol][iRow + 1] += pdFrom[M12];
1914  ppdColsm1[iCol][iRow + 2] += pdFrom[M13];
1915 
1916  ppdColsm1[iCol + 1][iRow] += pdFrom[M21];
1917  ppdColsm1[iCol + 1][iRow + 1] += pdFrom[M22];
1918  ppdColsm1[iCol + 1][iRow + 2] += pdFrom[M23];
1919 
1920  ppdColsm1[iCol + 2][iRow] += pdFrom[M31];
1921  ppdColsm1[iCol + 2][iRow + 1] += pdFrom[M32];
1922  ppdColsm1[iCol + 2][iRow + 2] += pdFrom[M33];
1923 }
1924 
1925 
1926 /* sottrae una matrice di tipo Mat3x3 in una data posizione
1927  * analoga ala precedente, con il meno (per evitare temporanei ecc) */
1928 void
1929 FullMatrixHandler::Sub(integer iRow, integer iCol, const Mat3x3& m)
1930 {
1931  /* iRow e iCol sono gli indici effettivi di riga e colonna
1932  * es. per il primo coefficiente:
1933  * iRow = 1, iCol = 1
1934  */
1935 
1936 #ifdef DEBUG
1937  IsValid();
1938 
1939  ASSERT(iRow > 0);
1940  ASSERT(iRow <= iNumRows - 2);
1941  ASSERT(iCol > 0);
1942  ASSERT(iCol <= iNumCols - 2);
1943 #endif /* DEBUG */
1944 
1945  const doublereal* pdFrom = m.pGetMat();
1946 
1947  ppdColsm1[iCol][iRow] -= pdFrom[M11];
1948  ppdColsm1[iCol][iRow + 1] -= pdFrom[M21];
1949  ppdColsm1[iCol][iRow + 2] -= pdFrom[M31];
1950 
1951  ppdColsm1[iCol + 1][iRow] -= pdFrom[M12];
1952  ppdColsm1[iCol + 1][iRow + 1] -= pdFrom[M22];
1953  ppdColsm1[iCol + 1][iRow + 2] -= pdFrom[M32];
1954 
1955  ppdColsm1[iCol + 2][iRow] -= pdFrom[M13];
1956  ppdColsm1[iCol + 2][iRow + 1] -= pdFrom[M23];
1957  ppdColsm1[iCol + 2][iRow + 2] -= pdFrom[M33];
1958 
1959 #if 0
1960  for (unsigned c = 0; c < 2; c++) {
1961  for (unsigned r = 0; r < 2; r++) {
1962  ppdColsm1[iCol + c][iRow + r] -= pdFrom[r];
1963  }
1964  pdFrom += 3;
1965  }
1966 #endif
1967 }
1968 
1969 void
1970 FullMatrixHandler::SubT(integer iRow, integer iCol, const Mat3x3& m)
1971 {
1972  /* iRow e iCol sono gli indici effettivi di riga e colonna
1973  * es. per il primo coefficiente:
1974  * iRow = 1, iCol = 1
1975  */
1976 
1977 #ifdef DEBUG
1978  IsValid();
1979 
1980  ASSERT(iRow > 0);
1981  ASSERT(iRow <= iNumRows - 2);
1982  ASSERT(iCol > 0);
1983  ASSERT(iCol <= iNumCols - 2);
1984 #endif /* DEBUG */
1985 
1986  const doublereal* pdFrom = m.pGetMat();
1987 
1988  ppdColsm1[iCol][iRow] -= pdFrom[M11];
1989  ppdColsm1[iCol][iRow + 1] -= pdFrom[M12];
1990  ppdColsm1[iCol][iRow + 2] -= pdFrom[M13];
1991 
1992  ppdColsm1[iCol + 1][iRow] -= pdFrom[M21];
1993  ppdColsm1[iCol + 1][iRow + 1] -= pdFrom[M22];
1994  ppdColsm1[iCol + 1][iRow + 2] -= pdFrom[M23];
1995 
1996  ppdColsm1[iCol + 2][iRow] -= pdFrom[M31];
1997  ppdColsm1[iCol + 2][iRow + 1] -= pdFrom[M32];
1998  ppdColsm1[iCol + 2][iRow + 2] -= pdFrom[M33];
1999 }
2000 
2001 
2002 /* mette una matrice di tipo Mat3x3 in una data posizione;
2003  * analoga alle precedenti */
2004 void
2005 FullMatrixHandler::Put(integer iRow, integer iCol, const Mat3x3& m)
2006 {
2007  /* iRow e iCol sono gli indici effettivi di riga e colonna
2008  * es. per il primo coefficiente:
2009  * iRow = 1, iCol = 1
2010  */
2011 
2012 #ifdef DEBUG
2013  IsValid();
2014 
2015  ASSERT(iRow > 0);
2016  ASSERT(iRow <= iNumRows - 2);
2017  ASSERT(iCol > 0);
2018  ASSERT(iCol <= iNumCols - 2);
2019 #endif /* DEBUG */
2020 
2021  const doublereal* pdFrom = m.pGetMat();
2022 
2023  ppdColsm1[iCol][iRow] = pdFrom[M11];
2024  ppdColsm1[iCol][iRow + 1] = pdFrom[M21];
2025  ppdColsm1[iCol][iRow + 2] = pdFrom[M31];
2026 
2027  ppdColsm1[iCol + 1][iRow] = pdFrom[M12];
2028  ppdColsm1[iCol + 1][iRow + 1] = pdFrom[M22];
2029  ppdColsm1[iCol + 1][iRow + 2] = pdFrom[M32];
2030 
2031  ppdColsm1[iCol + 2][iRow] = pdFrom[M13];
2032  ppdColsm1[iCol + 2][iRow + 1] = pdFrom[M23];
2033  ppdColsm1[iCol + 2][iRow + 2] = pdFrom[M33];
2034 
2035 #if 0
2036  for (unsigned c = 0; c < 2; c++) {
2037  for (unsigned r = 0; r < 2; r++) {
2038  ppdColsm1[iCol + c][iRow + r] = pdFrom[r];
2039  }
2040  pdFrom += 3;
2041  }
2042 #endif
2043 }
2044 
2045 void
2046 FullMatrixHandler::PutT(integer iRow, integer iCol, const Mat3x3& m)
2047 {
2048  /* iRow e iCol sono gli indici effettivi di riga e colonna
2049  * es. per il primo coefficiente:
2050  * iRow = 1, iCol = 1
2051  */
2052 
2053 #ifdef DEBUG
2054  IsValid();
2055 
2056  ASSERT(iRow > 0);
2057  ASSERT(iRow <= iNumRows - 2);
2058  ASSERT(iCol > 0);
2059  ASSERT(iCol <= iNumCols - 2);
2060 #endif /* DEBUG */
2061 
2062  const doublereal* pdFrom = m.pGetMat();
2063 
2064  ppdColsm1[iCol][iRow] = pdFrom[M11];
2065  ppdColsm1[iCol][iRow + 1] = pdFrom[M12];
2066  ppdColsm1[iCol][iRow + 2] = pdFrom[M13];
2067 
2068  ppdColsm1[iCol + 1][iRow] = pdFrom[M21];
2069  ppdColsm1[iCol + 1][iRow + 1] = pdFrom[M22];
2070  ppdColsm1[iCol + 1][iRow + 2] = pdFrom[M23];
2071 
2072  ppdColsm1[iCol + 2][iRow] = pdFrom[M31];
2073  ppdColsm1[iCol + 2][iRow + 1] = pdFrom[M32];
2074  ppdColsm1[iCol + 2][iRow + 2] = pdFrom[M33];
2075 }
2076 
2077 
2078 /* somma una matrice di tipo Mat3xN in una data posizione */
2079 void
2080 FullMatrixHandler::Add(integer iRow, integer iCol, const Mat3xN& m)
2081 {
2082  /* iRow e iCol sono gli indici effettivi di riga e colonna
2083  * es. per il primo coefficiente:
2084  * iRow = 1, iCol = 1 */
2085 
2086 #ifdef DEBUG
2087  IsValid();
2088 
2089  ASSERT(iRow > 0);
2090  ASSERT(iRow <= iNumRows - 2);
2091  ASSERT(iCol > 0);
2092  ASSERT(iCol <= iNumCols - m.iGetNumCols() + 1);
2093 #endif /* DEBUG */
2094 
2095  --iRow;
2096  --iCol;
2097  for (int r = 3; r > 0; r--) {
2098  for (integer c = m.iGetNumCols(); c > 0; c--) {
2099  ppdColsm1[iCol + c][iRow + r] += m(r, c);
2100  }
2101  }
2102 }
2103 
2104 /* sottrae una matrice di tipo Mat3xN in una data posizione */
2105 void
2106 FullMatrixHandler::Sub(integer iRow, integer iCol, const Mat3xN& m)
2107 {
2108  /* iRow e iCol sono gli indici effettivi di riga e colonna
2109  * es. per il primo coefficiente:
2110  * iRow = 1, iCol = 1 */
2111 
2112 #ifdef DEBUG
2113  IsValid();
2114 
2115  ASSERT(iRow > 0);
2116  ASSERT(iRow <= iNumRows - 2);
2117  ASSERT(iCol > 0);
2118  ASSERT(iCol <= iNumCols - m.iGetNumCols() + 1);
2119 #endif /* DEBUG */
2120 
2121  --iRow;
2122  --iCol;
2123  for (int r = 3; r > 0; r--) {
2124  for (integer c = m.iGetNumCols(); c > 0; c--) {
2125  ppdColsm1[iCol + c][iRow + r] -= m(r, c);
2126  }
2127  }
2128 }
2129 
2130 /* setta una matrice di tipo Mat3xN in una data posizione */
2131 void
2132 FullMatrixHandler::Put(integer iRow, integer iCol, const Mat3xN& m)
2133 {
2134  /* iRow e iCol sono gli indici effettivi di riga e colonna
2135  * es. per il primo coefficiente:
2136  * iRow = 1, iCol = 1 */
2137 
2138 #ifdef DEBUG
2139  IsValid();
2140 
2141  ASSERT(iRow > 0);
2142  ASSERT(iRow <= iNumRows - 2);
2143  ASSERT(iCol > 0);
2144  ASSERT(iCol <= iNumCols - m.iGetNumCols() + 1);
2145 #endif /* DEBUG */
2146 
2147  --iRow;
2148  --iCol;
2149  for (int r = 3; r > 0; r--) {
2150  for (integer c = m.iGetNumCols(); c > 0; c--) {
2151  ppdColsm1[iCol + c][iRow + r] = m(r, c);
2152  }
2153  }
2154 }
2155 
2156 /* somma una matrice di tipo Mat3xN in una data posizione, trasposta */
2157 void
2158 FullMatrixHandler::AddT(integer iRow, integer iCol, const Mat3xN& m)
2159 {
2160  /* iRow e iCol sono gli indici effettivi di riga e colonna
2161  * es. per il primo coefficiente:
2162  * iRow = 1, iCol = 1 */
2163 
2164 #ifdef DEBUG
2165  IsValid();
2166 
2167  ASSERT(iRow > 0);
2168  ASSERT(iRow <= iNumRows - m.iGetNumCols() + 1);
2169  ASSERT(iCol > 0);
2170  ASSERT(iCol <= iNumCols - 2);
2171 #endif /* DEBUG */
2172 
2173  --iRow;
2174  --iCol;
2175  for (int r = m.iGetNumCols(); r > 0; r--) {
2176  for (integer c = 3; c > 0; c--) {
2177  ppdColsm1[iCol + c][iRow + r] += m(c, r);
2178  }
2179  }
2180 }
2181 
2182 /* sottrae una matrice di tipo Mat3xN in una data posizione, trasposta */
2183 void
2184 FullMatrixHandler::SubT(integer iRow, integer iCol, const Mat3xN& m)
2185 {
2186  /* iRow e iCol sono gli indici effettivi di riga e colonna
2187  * es. per il primo coefficiente:
2188  * iRow = 1, iCol = 1 */
2189 
2190 #ifdef DEBUG
2191  IsValid();
2192 
2193  ASSERT(iRow > 0);
2194  ASSERT(iRow <= iNumRows - m.iGetNumCols() + 1);
2195  ASSERT(iCol > 0);
2196  ASSERT(iCol <= iNumCols - 2);
2197 #endif /* DEBUG */
2198 
2199  --iRow;
2200  --iCol;
2201  for (int r = m.iGetNumCols(); r > 0; r--) {
2202  for (integer c = 3; c > 0; c--) {
2203  ppdColsm1[iCol + c][iRow + r] -= m(c, r);
2204  }
2205  }
2206 }
2207 
2208 /* setta una matrice di tipo Mat3xN in una data posizione, trasposta */
2209 void
2210 FullMatrixHandler::PutT(integer iRow, integer iCol, const Mat3xN& m)
2211 {
2212  /* iRow e iCol sono gli indici effettivi di riga e colonna
2213  * es. per il primo coefficiente:
2214  * iRow = 1, iCol = 1 */
2215 
2216 #ifdef DEBUG
2217  IsValid();
2218 
2219  ASSERT(iRow > 0);
2220  ASSERT(iRow <= iNumRows - m.iGetNumCols() + 1);
2221  ASSERT(iCol > 0);
2222  ASSERT(iCol <= iNumCols - 2);
2223 #endif /* DEBUG */
2224 
2225  --iRow;
2226  --iCol;
2227  for (int r = m.iGetNumCols(); r > 0; r--) {
2228  for (integer c = 3; c > 0; c--) {
2229  ppdColsm1[iCol + c][iRow + r] = m(c, r);
2230  }
2231  }
2232 }
2233 
2234 /* somma una matrice di tipo MatNx3 in una data posizione */
2235 void
2236 FullMatrixHandler::Add(integer iRow, integer iCol, const MatNx3& m)
2237 {
2238 #ifdef DEBUG
2239  IsValid();
2240 
2241  ASSERT(iRow > 0);
2242  ASSERT(iRow <= iNumRows - m.iGetNumRows() + 1);
2243  ASSERT(iCol > 0);
2244  ASSERT(iCol <= iNumCols - 2);
2245 #endif /* DEBUG */
2246 
2247  --iRow;
2248  --iCol;
2249  for (int r = m.iGetNumRows(); r > 0; r--) {
2250  for (integer c = 3; c > 0; c--) {
2251  ppdColsm1[iCol + c][iRow + r] += m(r, c);
2252  }
2253  }
2254 }
2255 
2256 /* sottrae una matrice di tipo MatNx3 in una data posizione */
2257 void
2258 FullMatrixHandler::Sub(integer iRow, integer iCol, const MatNx3& m)
2259 {
2260 #ifdef DEBUG
2261  IsValid();
2262 
2263  ASSERT(iRow > 0);
2264  ASSERT(iRow <= iNumRows - m.iGetNumRows() + 1);
2265  ASSERT(iCol > 0);
2266  ASSERT(iCol <= iNumCols - 2);
2267 #endif /* DEBUG */
2268 
2269  --iRow;
2270  --iCol;
2271  for (int r = m.iGetNumRows(); r > 0; r--) {
2272  for (integer c = 3; c > 0; c--) {
2273  ppdColsm1[iCol + c][iRow + r] -= m(r, c);
2274  }
2275  }
2276 }
2277 
2278 /* setta una matrice di tipo MatNx3 in una data posizione */
2279 void
2280 FullMatrixHandler::Put(integer iRow, integer iCol, const MatNx3& m)
2281 {
2282 #ifdef DEBUG
2283  IsValid();
2284 
2285  ASSERT(iRow > 0);
2286  ASSERT(iRow <= iNumRows - m.iGetNumRows() + 1);
2287  ASSERT(iCol > 0);
2288  ASSERT(iCol <= iNumCols - 2);
2289 #endif /* DEBUG */
2290 
2291  --iRow;
2292  --iCol;
2293  for (int r = m.iGetNumRows(); r > 0; r--) {
2294  for (integer c = 3; c > 0; c--) {
2295  ppdColsm1[iCol + c][iRow + r] = m(r, c);
2296  }
2297  }
2298 }
2299 
2300 void
2301 FullMatrixHandler::CopyMatrixRow(integer dest_row,
2302  const FullMatrixHandler & source, integer source_row)
2303 {
2304  integer nc = iGetNumCols();
2305 
2306  ASSERT(dest_row >= 1);
2307  ASSERT(dest_row <= iGetNumRows());
2308  ASSERT(source_row >= 1);
2309  ASSERT(source_row <= source.iGetNumRows());
2310  ASSERT(nc == source.iGetNumCols());
2311 
2312  for (integer ic = 1; ic <= nc; ic++) {
2313  // dest(dest_row, ic) = source(source_row, ic);
2314  ppdColsm1[ic][dest_row] = source.ppdColsm1[ic][source_row];
2315  }
2316 }
2317 
2318 void
2319 FullMatrixHandler::CopyMatrixBlock(integer dest_row, integer dest_col,
2320  const FullMatrixHandler & source,
2321  integer source_start_row, integer source_end_row,
2322  integer source_start_col, integer source_end_col)
2323 {
2324  ASSERT(source_start_row >= 1);
2325  ASSERT(source_end_row <= source.iGetNumRows());
2326  ASSERT(source_start_row <= source_end_row);
2327  ASSERT(dest_row >= 1);
2328  ASSERT(dest_row + (source_end_row - source_start_row) <= iGetNumRows());
2329 
2330  ASSERT(source_start_col >= 1);
2331  ASSERT(source_end_col <= source.iGetNumCols());
2332  ASSERT(source_start_col <= source_end_col);
2333  ASSERT(dest_col >= 1);
2334  ASSERT(dest_col + (source_end_col - source_start_col) <= iGetNumCols());
2335 
2336  for (integer ir = source_start_row; ir <= source_end_row; ir++) {
2337  integer row = dest_row + (ir - source_start_row);
2338  for (integer ic = source_start_col; ic <= source_end_col; ic++) {
2339  integer col = dest_col + (ic - source_start_col);
2340  // dest(row, col) = source(ir, ic);
2341  ppdColsm1[col][row] = source.ppdColsm1[ic][ir];
2342  }
2343  }
2344 }
2345 
2346 void
2347 FullMatrixHandler::const_iterator::reset(bool is_end)
2348 {
2349  if (is_end) {
2350  elem.iRow = m.iNumRows;
2351  elem.iCol = m.iNumCols;
2352 
2353  } else {
2354  i_idx = 0;
2355  elem.iRow = 0;
2356  elem.iCol = 0;
2357  elem.dCoef = m.pdRaw[i_idx];
2358  }
2359 }
2360 
2361 FullMatrixHandler::const_iterator::const_iterator(const FullMatrixHandler& m, bool is_end)
2362 : m(m)
2363 {
2364  reset(is_end);
2365 }
2366 
2367 FullMatrixHandler::const_iterator::~const_iterator(void)
2368 {
2369  NO_OP;
2370 }
2371 
2372 const FullMatrixHandler::const_iterator&
2373 FullMatrixHandler::const_iterator::operator ++ (void) const
2374 {
2375 #if 0
2376  // NOTE: this version only iterates on non-zero entries
2377  do {
2378  ++i_idx;
2379  if (++elem.iRow == m.iNumRows) {
2380  if (++elem.iCol == m.iNumCols) {
2381  return *this;
2382  }
2383  elem.iRow = 0;
2384  }
2385  } while (m.pdRaw[i_idx] == 0.);
2386 #else
2387  // NOTE: this version iterates on all coefficients
2388  ++i_idx;
2389  if (++elem.iRow == m.iNumRows) {
2390  if (++elem.iCol == m.iNumCols) {
2391  return *this;
2392  }
2393  elem.iRow = 0;
2394  }
2395 #endif
2396 
2397  elem.dCoef = m.pdRaw[i_idx];
2398 
2399  return *this;
2400 }
2401 
2402 const SparseMatrixHandler::SparseMatrixElement *
2403 FullMatrixHandler::const_iterator::operator -> (void) const
2404 {
2405  return &elem;
2406 }
2407 
2408 const SparseMatrixHandler::SparseMatrixElement&
2409 FullMatrixHandler::const_iterator::operator * (void) const
2410 {
2411  return elem;
2412 }
2413 
2414 bool
2415 FullMatrixHandler::const_iterator::operator == (const FullMatrixHandler::const_iterator& op) const
2416 {
2417  return elem == op.elem;
2418 }
2419 
2420 bool
2421 FullMatrixHandler::const_iterator::operator != (const FullMatrixHandler::const_iterator& op) const
2422 {
2423  return elem != op.elem;
2424 }
2425 
2426 /* FullMatrixHandler - end */
2427 
FullMatrixHandler::PutT
void PutT(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1338
FullMatrixHandler::pdRaw
doublereal * pdRaw
Definition: fullmh.h:69
FullMatrixHandler::MatMul
void MatMul(const FullMatrixHandler &m1, const FullMatrixHandler &m2)
Definition: fullmh.cc:433
FullMatrixHandler::const_iterator::reset
void reset(bool is_end=false)
Definition: fullmh.cc:2347
FullMatrixHandler::operator-=
virtual MatrixHandler & operator-=(const SubMatrixHandler &SubMH)
Definition: fullmh.cc:412
FullMatrixHandler::const_iterator::operator++
const FullMatrixHandler::const_iterator & operator++(void) const
Definition: fullmh.cc:2373
MatrixHandler::iGetNumCols
virtual integer iGetNumCols(void) const =0
MatrixHandler::IncCoef
virtual void IncCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: mh.cc:374
FullMatrixHandler::const_iterator::elem
SparseMatrixHandler::SparseMatrixElement elem
Definition: fullmh.h:83
MBDYN_EXCEPT_ARGS
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Vec3
Definition: matvec3.h:98
FullMatrixHandler::pdRawm1
doublereal * pdRawm1
Definition: fullmh.h:70
FullMatrixHandler::const_iterator
Definition: fullmh.h:77
FullMatrixHandler::const_iterator::operator==
bool operator==(const FullMatrixHandler::const_iterator &op) const
Definition: fullmh.cc:2415
FullMatrixHandler::iMaxCols
integer iMaxCols
Definition: fullmh.h:67
FullMatrixHandler::Put
void Put(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1113
M22
Definition: matvec3.h:59
FullMatrixHandler::const_iterator::operator->
const SparseMatrixHandler::SparseMatrixElement * operator->(void) const
Definition: fullmh.cc:2403
Write
std::ostream & Write(std::ostream &out, const FullMatrixHandler &m, const char *s, const char *s2)
Definition: fullmh.cc:376
SubMatrixHandler::AddTo
virtual MatrixHandler & AddTo(MatrixHandler &MH) const =0
FullMatrixHandler::Detach
void Detach(void)
Definition: fullmh.cc:268
FullMatrixHandler::const_iterator::i_idx
integer i_idx
Definition: fullmh.h:82
SAFEDELETEARR
#define SAFEDELETEARR(pnt)
Definition: mynewmem.h:713
FullMatrixHandler::const_iterator::const_iterator
const_iterator(const FullMatrixHandler &m, bool is_end=false)
Definition: fullmh.cc:2361
VectorHandler::IncCoef
virtual void IncCoef(integer iRow, const doublereal &dCoef)=0
FullMatrixHandler::iNumRows
integer iNumRows
Definition: fullmh.h:63
MatrixHandler::DecCoef
virtual void DecCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: mh.cc:379
FullMatrixHandler::CopyMatrixRow
void CopyMatrixRow(integer dest_row, const FullMatrixHandler &source, integer source_row)
Definition: fullmh.cc:2301
FullMatrixHandler::Add
void Add(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1039
MatNx3::iGetNumRows
integer iGetNumRows(void) const
Definition: matvec3n.h:431
NO_OP
#define NO_OP
Definition: myassert.h:74
V2
Definition: matvec3.h:50
FullMatrixHandler::operator+=
virtual MatrixHandler & operator+=(const SubMatrixHandler &SubMH)
Definition: fullmh.cc:405
FullMatrixHandler::Resize
virtual void Resize(integer iNewRows, integer iNewCols)
Definition: fullmh.cc:174
FullMatrixHandler::ppdCols
doublereal ** ppdCols
Definition: fullmh.h:71
VectorHandler::iGetSize
virtual integer iGetSize(void) const =0
MyVectorHandler::pdVecm1
doublereal * pdVecm1
Definition: vh.h:158
FullMatrixHandler::const_iterator::~const_iterator
~const_iterator(void)
Definition: fullmh.cc:2367
M21
Definition: matvec3.h:58
SparseMatrixHandler::SparseMatrixElement_base::iCol
integer iCol
Definition: spmh.h:59
V3
Definition: matvec3.h:51
M11
Definition: matvec3.h:55
SparseMatrixHandler::SparseMatrixElement_base::iRow
integer iRow
Definition: spmh.h:58
VectorHandler::DecCoef
virtual void DecCoef(integer iRow, const doublereal &dCoef)=0
FullMatrixHandler::~FullMatrixHandler
virtual ~FullMatrixHandler(void)
Definition: fullmh.cc:167
MatrixHandler::MatTMatMul_base
virtual MatrixHandler & MatTMatMul_base(void(MatrixHandler::*op)(integer iRow, integer iCol, const doublereal &dCoef), MatrixHandler &out, const MatrixHandler &in) const
Definition: mh.cc:183
SubMatrixHandler::SubFrom
virtual MatrixHandler & SubFrom(MatrixHandler &MH) const =0
M12
Definition: matvec3.h:56
FullMatrixHandler::const_iterator::operator!=
bool operator!=(const FullMatrixHandler::const_iterator &op) const
Definition: fullmh.cc:2421
M33
Definition: matvec3.h:63
M32
Definition: matvec3.h:62
M31
Definition: matvec3.h:61
false
Definition: mbdyn.h:76
operator<<
std::ostream & operator<<(std::ostream &out, const FullMatrixHandler &m)
Definition: fullmh.cc:352
VariableSubMatrixHandler::AddTo
MatrixHandler & AddTo(MatrixHandler &MH) const
Definition: submat.h:1261
FullMatrixHandler::iRawSize
integer iRawSize
Definition: fullmh.h:66
M13
Definition: matvec3.h:57
true
Definition: mbdyn.h:77
MatrixHandler::PutCoef
virtual void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: mh.cc:384
FullMatrixHandler::MatTVecMul_base
virtual VectorHandler & MatTVecMul_base(void(VectorHandler::*op)(integer iRow, const doublereal &dCoef), VectorHandler &out, const VectorHandler &in) const
Definition: fullmh.cc:963
Mat3xN::iGetNumCols
integer iGetNumCols(void) const
Definition: matvec3n.h:298
FullMatrixHandler::iNumCols
integer iNumCols
Definition: fullmh.h:64
ASSERT
#define ASSERT(expression)
Definition: colamd.c:977
FullMatrixHandler::CopyMatrixBlock
void CopyMatrixBlock(integer dest_row, integer dest_col, const FullMatrixHandler &source, integer source_start_row, integer source_end_row, integer source_start_col, integer source_end_col)
Definition: fullmh.cc:2319
MatrixHandler::MatMatMul_base
virtual MatrixHandler & MatMatMul_base(void(MatrixHandler::*op)(integer iRow, integer iCol, const doublereal &dCoef), MatrixHandler &out, const MatrixHandler &in) const
Definition: mh.cc:135
FullMatrixHandler::CreateColRow
void CreateColRow(integer iNR, integer iNC)
Definition: fullmh.cc:58
FullMatrixHandler::MatTMatMul_base
MatrixHandler & MatTMatMul_base(void(MatrixHandler::*op)(integer iRow, integer iCol, const doublereal &dCoef), MatrixHandler &out, const MatrixHandler &in) const
Definition: fullmh.cc:744
VectorHandler::PutCoef
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
VariableSubMatrixHandler::SubFrom
MatrixHandler & SubFrom(MatrixHandler &MH) const
Definition: submat.h:1325
FullMatrixHandler::Sub
void Sub(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1076
defaultMemoryManager
#define defaultMemoryManager
Definition: mynewmem.h:691
FullMatrixHandler::iGetNumCols
virtual integer iGetNumCols(void) const
Definition: fullmh.h:229
FullMatrixHandler::const_iterator::operator*
const SparseMatrixHandler::SparseMatrixElement & operator*(void) const
Definition: fullmh.cc:2409
c
static std::stack< cleanup * > c
Definition: cleanup.cc:59
FullMatrixHandler::FullMatrixHandler
FullMatrixHandler(void)
Definition: fullmh.cc:119
FullMatrixHandler::MatMatMul_base
MatrixHandler & MatMatMul_base(void(MatrixHandler::*op)(integer iRow, integer iCol, const doublereal &dCoef), MatrixHandler &out, const MatrixHandler &in) const
Definition: fullmh.cc:601
Mat3x3::pGetMat
const doublereal * pGetMat(void) const
Definition: matvec3.h:743
FullMatrixHandler::Attach
void Attach(integer iNewRows, integer iNewCols, doublereal *pd, doublereal **ppd, integer iMSize=0, integer iMaxC=0)
Definition: fullmh.cc:290
V1
Definition: matvec3.h:49
Vec3::pGetVec
const doublereal * pGetVec(void) const
Definition: matvec3.h:192
M23
Definition: matvec3.h:60
FullMatrixHandler::const_iterator::m
const FullMatrixHandler & m
Definition: fullmh.h:81
SAFENEWARR
#define SAFENEWARR(pnt, item, sz)
Definition: mynewmem.h:701
FullMatrixHandler::bOwnsMemory
bool bOwnsMemory
Definition: fullmh.h:61
FullMatrixHandler::ppdColsm1
doublereal ** ppdColsm1
Definition: fullmh.h:72
FullMatrixHandler::SubT
void SubT(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1301
doublereal
double doublereal
Definition: colamd.c:52
integer
long int integer
Definition: colamd.c:51
FullMatrixHandler::AddT
void AddT(integer iRow, integer iCol, const FullMatrixHandler &source)
Definition: fullmh.cc:1264
FullMatrixHandler::MatVecMul_base
virtual VectorHandler & MatVecMul_base(void(VectorHandler::*op)(integer iRow, const doublereal &dCoef), VectorHandler &out, const VectorHandler &in) const
Definition: fullmh.cc:887
FullMatrixHandler::iGetNumRows
virtual integer iGetNumRows(void) const
Definition: fullmh.h:225
MatrixHandler::iGetNumRows
virtual integer iGetNumRows(void) const =0
FullMatrixHandler::Reset
void Reset(void)
Definition: fullmh.cc:44
FullMatrixHandler::operator=
FullMatrixHandler & operator=(const FullMatrixHandler &)
Definition: fullmh.cc:146