1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/eigjdqz.cc,v 1.9 2017/01/12 14:46:09 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 /*
35
36 3.1 User supplied subroutines
37 The user has to supply three problem dependent routines: one for the mul-
38 tiplication of a vector with the operator A, one for multiplication with B,
39 and one for performing the preconditioning operation. The subroutine to
40 multiply with A must be called AMUL and must have the following header:
41
42 subroutine AMUL( n, q, r )
43 c...............................................
44 c... Subroutine to compute r = Aq
45 c...............................................
46 integer n
47 double complex q(n), r(n)
48
49 q is the input vector, r the output vector. n is the dimension of the problem.
50 The subroutine to multiply with B must be called BMUL and must have the
51 following header:
52
53 subroutine BMUL( n, q, r )
54 c...............................................
55 c... Subroutine to compute r = Bq
56 c...............................................
57 integer n
58 double complex q(n), r(n)
59
60 Finally, the routine to perform the preconditioning operation must be called
61 PRECON and must have the header
62
63 subroutine PRECON( n, q )
64 c...............................................
65 c... Subroutine to compute q = K^-1 q
66 c...............................................
67 integer n
68 double complex q(n)
69
70 The preconditioning matrix should be an approximation of the matrix A -
71 B, with the prechosen target value. Preconditioning within the JDQZ
72 algorithm is described in section 3.4 of [1]. Preconditioning is not essential
73 for the correct behavior of the algorithm. It should improve the rate of
74 convergence, but leaving the vector q untouched should have no influence on
75 the correctness of the results.
76
77 */
78
79 #include "eigjdqz.h"
80
81 MBJDQZ* mbjdqzp;
82
83 extern "C" int
__FC_DECL__(amul)84 __FC_DECL__(amul)(integer *n, doublecomplex *q, doublecomplex *r)
85 {
86 ASSERT(mbjdqzp != 0);
87
88 mbjdqzp->AMul(*n, q, r);
89
90 return 0;
91 }
92
93 extern "C" int
__FC_DECL__(bmul)94 __FC_DECL__(bmul)(integer *n, doublecomplex *q, doublecomplex *r)
95 {
96 ASSERT(mbjdqzp != 0);
97
98 mbjdqzp->BMul(*n, q, r);
99
100 return 0;
101 }
102
103 extern "C" int
__FC_DECL__(precon)104 __FC_DECL__(precon)(integer *n, doublecomplex *q)
105 {
106 // leave untouched
107 return 0;
108 }
109
MBJDQZ(const NaiveMatrixHandler & A,const NaiveMatrixHandler & B)110 MBJDQZ::MBJDQZ(const NaiveMatrixHandler& A, const NaiveMatrixHandler& B)
111 : A(A), B(B), cnt(0)
112 {
113 NO_OP;
114 }
115
~MBJDQZ(void)116 MBJDQZ::~MBJDQZ(void)
117 {
118 NO_OP;
119 }
120
121 void
Mul(integer n,doublecomplex * q,doublecomplex * r,const NaiveMatrixHandler & M)122 MBJDQZ::Mul(integer n, doublecomplex *q, doublecomplex *r, const NaiveMatrixHandler& M)
123 {
124 ASSERT(n == M.iGetNumRows());
125 ASSERT(n == M.iGetNumCols());
126
127 for (integer i = 0; i < n; i++) {
128 r[i].r = 0.;
129 r[i].i = 0.;
130 }
131
132 for (NaiveMatrixHandler::const_iterator i = M.begin(); i != M.end(); ++i) {
133 r[i->iRow].r += i->dCoef*q[i->iCol].r;
134 r[i->iRow].i += i->dCoef*q[i->iCol].i;
135 }
136 }
137
138 void
AMul(integer n,doublecomplex * q,doublecomplex * r)139 MBJDQZ::AMul(integer n, doublecomplex *q, doublecomplex *r)
140 {
141 #define CNT (100)
142 if (!(cnt % CNT)) {
143 silent_cerr("\r" "cnt=" << cnt);
144 }
145
146 cnt++;
147
148 Mul(n, q, r, A);
149 }
150
151 void
BMul(integer n,doublecomplex * q,doublecomplex * r)152 MBJDQZ::BMul(integer n, doublecomplex *q, doublecomplex *r)
153 {
154 Mul(n, q, r, B);
155 }
156
157 unsigned
Cnt(void) const158 MBJDQZ::Cnt(void) const
159 {
160 return cnt;
161 }
162