1 /*
2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
3 * Applied Mathematics, Norway.
4 *
5 * Contact information: E-mail: tor.dokken@sintef.no
6 * SINTEF ICT, Department of Applied Mathematics,
7 * P.O. Box 124 Blindern,
8 * 0314 Oslo, Norway.
9 *
10 * This file is part of SISL.
11 *
12 * SISL is free software: you can redistribute it and/or modify
13 * it under the terms of the GNU Affero General Public License as
14 * published by the Free Software Foundation, either version 3 of the
15 * License, or (at your option) any later version.
16 *
17 * SISL is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Affero General Public License for more details.
21 *
22 * You should have received a copy of the GNU Affero General Public
23 * License along with SISL. If not, see
24 * <http://www.gnu.org/licenses/>.
25 *
26 * In accordance with Section 7(b) of the GNU Affero General Public
27 * License, a covered work must retain the producer line in every data
28 * file that is created or manipulated using SISL.
29 *
30 * Other Usage
31 * You can be released from the requirements of the license by purchasing
32 * a commercial license. Buying such a license is mandatory as soon as you
33 * develop commercial activities involving the SISL library without
34 * disclosing the source code of your own applications.
35 *
36 * This file may be used in accordance with the terms contained in a
37 * written agreement between you and SINTEF ICT.
38 */
39
40 #include "sisl-copyright.h"
41
42 /*
43 *
44 * $Id: s6mulvec.c,v 1.1 1994-04-21 12:10:42 boh Exp $
45 *
46 */
47
48
49 #define S6MULVEC
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s6mulvec(double ematrix[],double evect[],double eright[])55 s6mulvec (double ematrix[], double evect[], double eright[])
56 #else
57 void
58 s6mulvec (ematrix, evect, eright)
59 double ematrix[];
60 double evect[];
61 double eright[];
62 #endif
63 /*
64 *************************************************************************
65 *
66 * Purpose: To multiply a 4*4 matrix by a 3-D vector.
67 *
68 * Input:
69 * Ematrix - The matrix (length 16).
70 * Evect - The vector (length 3).
71 *
72 * Output:
73 * Eright - The resulting 3-D vector.
74 *
75 * Written by: A.M. Ytrehus, SI Oslo Nov.91.
76 * After FORTRAN (P6mvec), written by: S. Meen SI.
77 *****************************************************************
78 */
79 {
80 double svect[4];
81 double sright[4];
82 double tdum;
83 int ki;
84 int kj;
85
86 /* Store the 3D evect-array in the 4D svect-array. */
87
88 for (ki = 0; ki < 3; ki++)
89 svect[ki] = evect[ki];
90
91 svect[3] = (double) 1.0;
92
93
94 /* Multiply the matrix by the vector svect.
95 The result is stored in the 4-D local vector sright.*/
96
97 for (ki = 0; ki < 4; ki++)
98 {
99 tdum = (double) 0.0;
100
101 for (kj = 0; kj < 4; kj++)
102 tdum += ematrix[4 * ki + kj] * svect[kj];
103
104 sright[ki] = tdum;
105 }
106
107
108 /*
109 * Check if the bottom line of the matrix is 0,0,0,1.
110 * In that case, just store the first three elements of svect in evect.
111 * Else, divide all 3 first elements of svect by svect[3].
112 * --------------------------------------------------------------------
113 */
114
115 for (ki = 0; ki < 3; ki++)
116 eright[ki] = sright[ki];
117
118 if (!(ematrix[12] == (double) 0.0) && (ematrix[13] == (double) 0.0) &&
119 (ematrix[14] == (double) 0.0) && (ematrix[15] == (double) 1.0))
120 {
121 tdum = sright[3];
122
123 for (ki = 0; ki < 3; ki++)
124 eright[ki] = sright[ki] / tdum;
125 }
126
127 goto out;
128
129 out:
130 return;
131 }
132