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: s6mvec.c,v 1.2 2001-03-19 15:59:02 afr Exp $
45 *
46 */
47
48
49 #define S6MVEC
50
51 #include "sislP.h"
52
53 #if defined(SISLNEEDPROTOTYPES)
54 void
s6mvec(double emat[],double evec1[],int inbvec,double evec2[])55 s6mvec(double emat[],double evec1[],int inbvec,double evec2[])
56 #else
57 void s6mvec(emat,evec1,inbvec,evec2)
58 double emat[];
59 double evec1[];
60 int inbvec;
61 double evec2[];
62 #endif
63 /*
64 *********************************************************************
65 *
66 * PURPOSE : To post multiply a matrix by inbvec 3-D vectors.
67 *
68 * INPUT : emat - The matrix (16 elements, Homogenous coordinates)
69 * evec1 - The input vectors (3-D coordinates)
70 * inbvec - The number of input vectors
71 *
72 * OUTPUT : evec2 - The resulting vectors (3-D coordinates)
73 *
74 *-
75 * CALLS :
76 *
77 * WRITTEN BY : Tor Dokken, SI, Oslo, Norway. 1988-may-24
78 *
79 *********************************************************************
80 */
81 {
82 double svec2[3]; /* Temporary vector */
83 double *svec1; /* Pointer to current vector */
84 double *svec3; /* pointer to current result vector */
85 register double tdum; /* Temporary storage of real */
86 register int ki,kj,kl,kp; /* Control variables in loops */
87 int kstop; /* Stop condition for loop */
88
89 /* Multiply matrix by all vectors */
90
91 kstop = 3*inbvec;
92
93 for (kl=0;kl<kstop;kl=kl+3)
94 {
95 /* Multiply rotational part of the matrix by the vector */
96
97 svec1 = evec1 + kl;
98 svec3 = evec2 + kl;
99
100 for (ki=0;ki<3;ki++)
101 {
102 kp = ki;
103
104 tdum = DZERO;
105 for (kj=0;kj<3;kj++)
106 {
107 tdum += emat[kp]*svec1[kj];
108 kp += 4;
109 }
110
111 /* Add translation part */
112 svec2[ki] = tdum + emat[kp];
113 }
114 /* Check if the bottom row is 0,0,0,1 */
115
116 if (DNEQUAL(emat[3],DZERO) || DNEQUAL(emat[7],DZERO) ||
117 DNEQUAL(emat[11],DZERO) || DNEQUAL(emat[15],(double)1.0))
118 {
119 /* Compute last element of vector */
120
121 tdum = evec1[0]*emat[3] + evec1[1]*emat[7] + evec1[2]*emat[11];
122 if (DNEQUAL(tdum,DZERO))
123 {
124 for (ki=0;ki<3;ki++)
125 svec2[ki] /= tdum;
126 }
127 }
128 svec3[0] = svec2[0];
129 svec3[1] = svec2[1];
130 svec3[2] = svec2[2];
131 }
132
133 return;
134 }
135