1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2008-2008 - INRIA - Antoine ELIAS <antoine.elias@scilab.org>
4 *
5 * Copyright (C) 2012 - 2016 - Scilab Enterprises
6 *
7 * This file is hereby licensed under the terms of the GNU GPL v2.0,
8 * pursuant to article 5.3.4 of the CeCILL v.2.1.
9 * This file was originally licensed under the terms of the CeCILL v2.1,
10 * and continues to be available under such terms.
11 * For more information, see the COPYING file which you should have received
12 * along with this program.
13 *
14 */
15
16 #include <string.h> // for memset fonction
17 #include "elem_common.h"
18 #include "matrix_multiplication.h"
19
iMultiComplexMatrixByComplexMatrix(double * _pdblReal1,double * _pdblImg1,int _iRows1,int _iCols1,double * _pdblReal2,double * _pdblImg2,int _iRows2,int _iCols2,double * _pdblRealOut,double * _pdblImgOut)20 int iMultiComplexMatrixByComplexMatrix(
21 double *_pdblReal1, double *_pdblImg1, int _iRows1, int _iCols1,
22 double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2,
23 double *_pdblRealOut, double *_pdblImgOut)
24 {
25 double dblOne = 1;
26 double dblMinusOne = -1;
27 double dblZero = 0;
28
29 int iRowsOut = _iRows1;
30
31 //Cr <- 1*Ar*Br + 0*Cr
32 F2C(dgemm)("N", "N", &iRowsOut, &_iCols2, &_iCols1, &dblOne , _pdblReal1 , &_iRows1, _pdblReal2 , &_iRows2, &dblZero , _pdblRealOut , &iRowsOut);
33 //Cr <- -1*Ai*Bi + 1*Cr
34 F2C(dgemm)("N", "N", &iRowsOut, &_iCols2, &_iCols1, &dblMinusOne , _pdblImg1 , &_iRows1, _pdblImg2 , &_iRows2, &dblOne , _pdblRealOut , &iRowsOut);
35 //Ci <- 1*Ar*Bi + 0*Ci
36 F2C(dgemm)("N", "N", &iRowsOut, &_iCols2, &_iCols1, &dblOne , _pdblReal1 , &_iRows1, _pdblImg2 , &_iRows2, &dblZero , _pdblImgOut , &iRowsOut);
37 //Ci <- 1*Ai*Br + 1*Ci
38 F2C(dgemm)("N", "N", &iRowsOut, &_iCols2, &_iCols1, &dblOne , _pdblImg1 , &_iRows1, _pdblReal2 , &_iRows2, &dblOne , _pdblImgOut , &iRowsOut);
39
40 return 0;
41 }
42
iMultiRealMatrixByRealMatrix(double * _pdblReal1,int _iRows1,int _iCols1,double * _pdblReal2,int _iRows2,int _iCols2,double * _pdblRealOut)43 int iMultiRealMatrixByRealMatrix(
44 double *_pdblReal1, int _iRows1, int _iCols1,
45 double *_pdblReal2, int _iRows2, int _iCols2,
46 double *_pdblRealOut)
47 {
48 double dblOne = 1;
49 double dblZero = 0;
50
51 C2F(dgemm)("n", "n", &_iRows1, &_iCols2, &_iCols1, &dblOne,
52 _pdblReal1 , &_iRows1 ,
53 _pdblReal2, &_iRows2, &dblZero,
54 _pdblRealOut , &_iRows1);
55 return 0;
56 }
57
iMultiRealMatrixByComplexMatrix(double * _pdblReal1,int _iRows1,int _iCols1,double * _pdblReal2,double * _pdblImg2,int _iRows2,int _iCols2,double * _pdblRealOut,double * _pdblImgOut)58 int iMultiRealMatrixByComplexMatrix(
59 double *_pdblReal1, int _iRows1, int _iCols1,
60 double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2,
61 double *_pdblRealOut, double *_pdblImgOut)
62 {
63 double dblOne = 1;
64 double dblZero = 0;
65 iMultiRealMatrixByRealMatrix(_pdblReal1, _iRows1, _iCols1, _pdblReal2, _iRows2, _iCols2, _pdblRealOut);
66 C2F(dgemm)("n", "n", &_iRows1, &_iCols2, &_iCols1, &dblOne,
67 _pdblReal1 , &_iRows1 ,
68 _pdblImg2, &_iRows2, &dblZero,
69 _pdblImgOut , &_iRows1);
70
71 return 0;
72 }
73
iMultiComplexMatrixByRealMatrix(double * _pdblReal1,double * _pdblImg1,int _iRows1,int _iCols1,double * _pdblReal2,int _iRows2,int _iCols2,double * _pdblRealOut,double * _pdblImgOut)74 int iMultiComplexMatrixByRealMatrix(
75 double *_pdblReal1, double *_pdblImg1, int _iRows1, int _iCols1,
76 double *_pdblReal2, int _iRows2, int _iCols2,
77 double *_pdblRealOut, double *_pdblImgOut)
78 {
79 double dblOne = 1;
80 double dblZero = 0;
81 iMultiRealMatrixByRealMatrix(_pdblReal1, _iRows1, _iCols1, _pdblReal2, _iRows2, _iCols2, _pdblRealOut);
82 C2F(dgemm)("n", "n", &_iRows1, &_iCols2, &_iCols1, &dblOne,
83 _pdblImg1 , &_iRows1 ,
84 _pdblReal2, &_iRows2, &dblZero,
85 _pdblImgOut , &_iRows1);
86 return 0;
87 }
88
89
iMultiRealScalarByRealMatrix(double _dblReal1,double * _pdblReal2,int _iRows2,int _iCols2,double * _pdblRealOut)90 int iMultiRealScalarByRealMatrix(
91 double _dblReal1,
92 double *_pdblReal2, int _iRows2, int _iCols2,
93 double *_pdblRealOut)
94 {
95 int iOne = 1;
96 int iSize2 = _iRows2 * _iCols2;
97
98 C2F(dcopy)(&iSize2, _pdblReal2, &iOne, _pdblRealOut, &iOne);
99 C2F(dscal)(&iSize2, &_dblReal1, _pdblRealOut, &iOne);
100 return 0;
101 }
102
iMultiRealScalarByComplexMatrix(double _dblReal1,double * _pdblReal2,double * _pdblImg2,int _iRows2,int _iCols2,double * _pdblRealOut,double * _pdblImgOut)103 int iMultiRealScalarByComplexMatrix(
104 double _dblReal1,
105 double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2,
106 double *_pdblRealOut, double *_pdblImgOut)
107 {
108 int iOne = 1;
109 int iSize2 = _iRows2 * _iCols2;
110
111 C2F(dcopy)(&iSize2, _pdblReal2, &iOne, _pdblRealOut, &iOne);
112 C2F(dcopy)(&iSize2, _pdblImg2, &iOne, _pdblImgOut, &iOne);
113 C2F(dscal)(&iSize2, &_dblReal1, _pdblRealOut, &iOne);
114 C2F(dscal)(&iSize2, &_dblReal1, _pdblImgOut, &iOne);
115 return 0;
116 }
117
iMultiComplexScalarByRealMatrix(double _dblReal1,double _dblImg1,double * _pdblReal2,int _iRows2,int _iCols2,double * _pdblRealOut,double * _pdblImgOut)118 int iMultiComplexScalarByRealMatrix(
119 double _dblReal1, double _dblImg1,
120 double *_pdblReal2, int _iRows2, int _iCols2,
121 double *_pdblRealOut, double *_pdblImgOut)
122 {
123 int iOne = 1;
124 int iSize2 = _iRows2 * _iCols2;
125
126 C2F(dcopy)(&iSize2, _pdblReal2, &iOne, _pdblRealOut, &iOne);
127 C2F(dcopy)(&iSize2, _pdblReal2, &iOne, _pdblImgOut, &iOne);
128 C2F(dscal)(&iSize2, &_dblReal1, _pdblRealOut, &iOne);
129 C2F(dscal)(&iSize2, &_dblImg1, _pdblImgOut, &iOne);
130 return 0;
131 }
132
iMultiComplexScalarByComplexMatrix(double _dblReal1,double _dblImg1,double * _pdblReal2,double * _pdblImg2,int _iRows2,int _iCols2,double * _pdblRealOut,double * _pdblImgOut)133 int iMultiComplexScalarByComplexMatrix(
134 double _dblReal1, double _dblImg1,
135 double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2,
136 double *_pdblRealOut, double *_pdblImgOut)
137 {
138 int iOne = 1;
139 int iSize2 = _iRows2 * _iCols2;
140 doublecomplex *pDataIn1 =
141 oGetDoubleComplexFromPointer(&_dblReal1, &_dblImg1, 1);
142 doublecomplex *pDataIn2 =
143 oGetDoubleComplexFromPointer(_pdblReal2, _pdblImg2, _iRows2 * _iCols2);
144
145 C2F(zscal)(&iSize2, pDataIn1, pDataIn2, &iOne);
146
147 vGetPointerFromDoubleComplex(pDataIn2, iSize2, _pdblRealOut, _pdblImgOut);
148 vFreeDoubleComplexFromPointer(pDataIn2);
149 vFreeDoubleComplexFromPointer(pDataIn1);
150 return 0;
151 }
152
153
154 /* (a1 + a2 * %s + ... ) * (a1 + a2 * %s + ... )*/
iMultiScilabPolynomByScilabPolynom(double * _pdblReal1,int _iRank1,double * _pdblReal2,int _iRank2,double * _pdblRealOut,int _iRankOut)155 int iMultiScilabPolynomByScilabPolynom(
156 double *_pdblReal1, int _iRank1,
157 double *_pdblReal2, int _iRank2,
158 double *_pdblRealOut, int _iRankOut)
159 {
160 int i1 = 0;
161 int i2 = 0;
162 double dblMult = 0.0;
163 double dblAdd = 0.0;
164
165 memset(_pdblRealOut, 0x00, _iRankOut * sizeof(double));
166
167 for (i1 = 0 ; i1 < _iRank1 ; i1++)
168 {
169 for (i2 = 0 ; i2 < _iRank2 ; i2++)
170 {
171 dblMult = _pdblReal1[i1] * _pdblReal2[i2];
172 dblAdd = _pdblRealOut[i1 + i2] + dblMult;
173 if (fabs(dblAdd) > 2 * nc_eps() * Max(fabs(_pdblRealOut[i1 + i2]), fabs(dblMult)))
174 {
175 _pdblRealOut[i1 + i2] = dblAdd;
176 }
177 else
178 {
179 _pdblRealOut[i1 + i2] = 0.0;
180 }
181 }
182 }
183 return 0;
184 }
185
186 /* ((a1 +ib1) + (a2+ib2) * %s + ... ) (a1 + a2 * %s + ... ) */
iMultiComplexPolyByScilabPolynom(double * _pdblReal1,double * _pdblImg1,int _iRank1,double * _pdblReal2,int _iRank2,double * _pdblRealOut,double * _pdblImgOut,int _iRankOut)187 int iMultiComplexPolyByScilabPolynom(
188 double *_pdblReal1, double *_pdblImg1, int _iRank1,
189 double *_pdblReal2, int _iRank2,
190 double *_pdblRealOut, double *_pdblImgOut, int _iRankOut)
191 {
192 int i1 = 0;
193 int i2 = 0;
194
195 memset(_pdblRealOut , 0x00, _iRankOut * sizeof(double));
196 memset(_pdblImgOut , 0x00, _iRankOut * sizeof(double));
197 for (i1 = 0 ; i1 < _iRank1 ; i1++)
198 {
199 for (i2 = 0 ; i2 < _iRank2 ; i2++)
200 {
201 _pdblRealOut[i1 + i2] += _pdblReal1[i1] * _pdblReal2[i2];
202 _pdblImgOut[i1 + i2] += _pdblImg1[i1] * _pdblReal2[i2];
203 }
204 }
205 return 0;
206 }
207
208 /* (a1 + a2 * %s + ... ) * ((a1 +ib1) + (a2+ib2) * %s + ... )*/
iMultiScilabPolynomByComplexPoly(double * _pdblReal1,int _iRank1,double * _pdblReal2,double * _pdblImg2,int _iRank2,double * _pdblRealOut,double * _pdblImgOut,int _iRankOut)209 int iMultiScilabPolynomByComplexPoly(
210 double *_pdblReal1, int _iRank1,
211 double *_pdblReal2, double *_pdblImg2, int _iRank2,
212 double *_pdblRealOut, double *_pdblImgOut, int _iRankOut)
213 {
214 return iMultiComplexPolyByScilabPolynom(
215 _pdblReal2, _pdblImg2, _iRank2,
216 _pdblReal1, _iRank1,
217 _pdblRealOut, _pdblImgOut, _iRankOut);
218
219 }
220
221 /* ((a1 +ib1) + (a2+ib2) * %s + ... ) * ((a1 +ib1) + (a2+ib2) * %s + ... ) */
iMultiComplexPolyByComplexPoly(double * _pdblReal1,double * _pdblImg1,int _iRank1,double * _pdblReal2,double * _pdblImg2,int _iRank2,double * _pdblRealOut,double * _pdblImgOut,int _iRankOut)222 int iMultiComplexPolyByComplexPoly(
223 double *_pdblReal1, double *_pdblImg1, int _iRank1,
224 double *_pdblReal2, double *_pdblImg2, int _iRank2,
225 double *_pdblRealOut, double *_pdblImgOut, int _iRankOut)
226 {
227 int i1 = 0;
228 int i2 = 0;
229
230 memset(_pdblRealOut , 0x00, _iRankOut * sizeof(double));
231 memset(_pdblImgOut , 0x00, _iRankOut * sizeof(double));
232 for (i1 = 0 ; i1 < _iRank1 ; i1++)
233 {
234 for (i2 = 0 ; i2 < _iRank2 ; i2++)
235 {
236 _pdblRealOut[i1 + i2] += _pdblReal1[i1] * _pdblReal2[i2] - _pdblImg1[i1] * _pdblImg2[i2];
237 _pdblImgOut[i1 + i2] += _pdblImg1[i1] * _pdblReal2[i2] + _pdblImg2[i2] * _pdblReal1[i1];
238 }
239 }
240 return 0;
241 }
242
iDotMultiplyRealMatrixByRealMatrix(double * _pdblReal1,double * _pdblReal2,double * _pdblRealOut,int _iRowsOut,int _iColsOut)243 int iDotMultiplyRealMatrixByRealMatrix(
244 double* _pdblReal1,
245 double* _pdblReal2,
246 double* _pdblRealOut, int _iRowsOut, int _iColsOut)
247 {
248 int i = 0;
249
250 for (i = 0 ; i < _iRowsOut * _iColsOut ; i++)
251 {
252 _pdblRealOut[i] = _pdblReal1[i] * _pdblReal2[i];
253 }
254 return 0;
255 }
256
iDotMultiplyRealMatrixByComplexMatrix(double * _pdblReal1,double * _pdblReal2,double * _pdblImg2,double * _pdblRealOut,double * _pdblImgOut,int _iRowsOut,int _iColsOut)257 int iDotMultiplyRealMatrixByComplexMatrix(
258 double* _pdblReal1,
259 double* _pdblReal2, double* _pdblImg2,
260 double* _pdblRealOut, double* _pdblImgOut, int _iRowsOut, int _iColsOut)
261 {
262 int i = 0;
263
264 for (i = 0 ; i < _iRowsOut * _iColsOut ; i++)
265 {
266 _pdblRealOut[i] = _pdblReal1[i] * _pdblReal2[i];
267 _pdblImgOut[i] = _pdblReal1[i] * _pdblImg2[i];
268 }
269 return 0;
270 }
271
iDotMultiplyComplexMatrixByRealMatrix(double * _pdblReal1,double * _pdblImg1,double * _pdblReal2,double * _pdblRealOut,double * _pdblImgOut,int _iRowsOut,int _iColsOut)272 int iDotMultiplyComplexMatrixByRealMatrix(
273 double* _pdblReal1, double* _pdblImg1,
274 double* _pdblReal2,
275 double* _pdblRealOut, double* _pdblImgOut, int _iRowsOut, int _iColsOut)
276 {
277 int i = 0;
278
279 for (i = 0 ; i < _iRowsOut * _iColsOut ; i++)
280 {
281 _pdblRealOut[i] = _pdblReal1[i] * _pdblReal2[i];
282 _pdblImgOut[i] = _pdblImg1[i] * _pdblReal2[i];
283 }
284 return 0;
285 }
286
iDotMultiplyComplexMatrixByComplexMatrix(double * _pdblReal1,double * _pdblImg1,double * _pdblReal2,double * _pdblImg2,double * _pdblRealOut,double * _pdblImgOut,int _iRowsOut,int _iColsOut)287 int iDotMultiplyComplexMatrixByComplexMatrix(
288 double* _pdblReal1, double* _pdblImg1,
289 double* _pdblReal2, double* _pdblImg2,
290 double* _pdblRealOut, double* _pdblImgOut, int _iRowsOut, int _iColsOut)
291 {
292 int i = 0;
293
294 for (i = 0 ; i < _iRowsOut * _iColsOut ; i++)
295 {
296 _pdblRealOut[i] = _pdblReal1[i] * _pdblReal2[i];
297 _pdblRealOut[i] -= _pdblImg1[i] * _pdblImg2[i];
298
299 _pdblImgOut[i] = _pdblImg1[i] * _pdblReal2[i];
300 _pdblImgOut[i] += _pdblReal1[i] * _pdblImg2[i];
301 }
302 return 0;
303 }
304
305