1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Antoine Elias
4 * Copyright (C) 2011 - DIGITEO - Cedric Delamarre
5 *
6 * Copyright (C) 2012 - 2016 - Scilab Enterprises
7 *
8 * This file is hereby licensed under the terms of the GNU GPL v2.0,
9 * pursuant to article 5.3.4 of the CeCILL v.2.1.
10 * This file was originally licensed under the terms of the CeCILL v2.1,
11 * and continues to be available under such terms.
12 * For more information, see the COPYING file which you should have received
13 * along with this program.
14 *
15 */
16
17 #include "matrix_kronecker.h"
18 #include "elem_common.h"
19
20 #include <stdio.h>
21 // Real .*. Real
vKronR(double * _pdblDataIn1,int _iIncIn1,int _iRowsIn1,int _iColsIn1,double * _pdblDataIn2,int _iIncIn2,int _iRowsIn2,int _iColsIn2,double * _pdblDataOut,int _iIncOut)22 void vKronR(double* _pdblDataIn1, int _iIncIn1, int _iRowsIn1, int _iColsIn1,
23 double* _pdblDataIn2, int _iIncIn2, int _iRowsIn2, int _iColsIn2,
24 double* _pdblDataOut, int _iIncOut)
25 {
26 //int iOne = 1;
27 //int iLoop1 = 0;
28 //int iLoop2 = 0;
29 //int iLoop3 = 0;
30 //int iIndex1 = - _iIncIn1;
31 //int iIndex2 = - _iColsIn2;
32
33 //for(iLoop1 = 0; iLoop1 < _iColsIn1; iLoop1++)
34 //{
35 // int iIndex3 = 0; //Index a 1 ???
36 // iIndex1 += _iIncIn1;
37 // iIndex2 += _iColsIn2;
38 // for(iLoop2 = 0; iLoop2 < _iColsIn2; iLoop2++)
39 // {
40 // int iIndex4 = iIndex1;
41 // int iIndex5 = (iLoop2 + iIndex2) * _iIncOut;
42 // for(iLoop3 = 0; iLoop3 < _iRowsIn1; iLoop3++)
43 // {
44 // C2F(dcopy)(&_iRowsIn2, &_pdblDataIn2[iIndex3], &iOne, &_pdblDataOut[iIndex5], &iOne);
45 // //ddscals(&_pdblDataOut[iIndex5], _iRowsIn2, _pdblDataIn1[iIndex4], &_pdblDataOut[iIndex5]);
46 // C2F(dscal)(&_iRowsIn2, &_pdblDataIn1[iIndex4], &_pdblDataOut[iIndex5], &iOne);
47
48 // iIndex5 += _iRowsIn2;
49 // iIndex4++;
50 // }
51 // iIndex3 += _iIncIn2;
52 // }
53 //}
54 int iSize1 = _iRowsIn1 * _iColsIn1;
55 int iSize2 = _iRowsIn2 * _iColsIn2;
56
57 int i = 0;
58 int j = 0;
59 int k = 0;
60 int l = 0;
61 int idx = 0;
62
63 for (i = 0 ; i < iSize1 ; i += _iRowsIn1)
64 {
65 int i2 = i + _iRowsIn1;
66 for (j = 0 ; j < iSize2 ; j += _iRowsIn2)
67 {
68 int j2 = j + _iRowsIn2;
69 for (k = i ; k < i2 ; k++)
70 {
71 //double c = _pdblDataIn1[k];
72 for (l = j ; l < j2 ; l++)
73 {
74 _pdblDataOut[idx] = _pdblDataIn1[k] * _pdblDataIn2[l];
75 idx++;
76 }
77 }
78 }
79 }
80 }
81
82 // Complex .*. Complex
vKronC(double * _pdblRealIn1,double * _pdblImgIn1,int _iIncIn1,int _iRowsIn1,int _iColsIn1,double * _pdblRealIn2,double * _pdblImgIn2,int _iIncIn2,int _iRowsIn2,int _iColsIn2,double * _pdblRealOut,double * _pdblImgOut,int _iIncOut)83 void vKronC(double* _pdblRealIn1, double* _pdblImgIn1, int _iIncIn1, int _iRowsIn1, int _iColsIn1,
84 double* _pdblRealIn2, double* _pdblImgIn2, int _iIncIn2, int _iRowsIn2, int _iColsIn2,
85 double* _pdblRealOut, double* _pdblImgOut, int _iIncOut)
86 {
87 //int iLoop1 = 0, iLoop2 = 0, iLoop3 = 0, iLoop4 = 0;
88 //int iIndex1 = -_iIncIn1;
89 //int iIndex2 = -_iColsIn2;
90 //for(iLoop1 = 0; iLoop1 < _iColsIn1; iLoop1++)
91 //{
92 // int iIndex3 = 0;
93 // iIndex1 += _iIncIn1;
94 // iIndex2 += _iColsIn2;
95 // for(iLoop2 = 0; iLoop2 < _iColsIn2; iLoop2++)
96 // {
97 // int iIndex4 = iIndex1;
98 // int iIndex5 = (iLoop2 + iIndex2) * _iIncOut;
99 // for(iLoop3 = 0; iLoop3 < _iRowsIn1; iLoop3++)
100 // {
101 // for(iLoop4 = 0; iLoop4 < _iRowsIn2; iLoop4++)
102 // {
103 // _pdblRealOut[iIndex5 + iLoop4] =
104 // _pdblRealIn1[iIndex4] * _pdblRealIn2[iIndex3 + iLoop4] -
105 // _pdblImgIn1[iIndex4] * _pdblImgIn2[iIndex3 + iLoop4];
106 // _pdblImgOut[iIndex5 + iLoop4] =
107 // _pdblRealIn1[iIndex4] * _pdblImgIn2[iIndex3 + iLoop4] +
108 // _pdblImgIn1[iIndex4] * _pdblRealIn2[iIndex3 + iLoop4];
109 // }
110 // iIndex5 += _iRowsIn2;
111 // iIndex4++;
112 // }
113 // iIndex3 += _iIncIn2;
114 // }
115 //}
116 int iSize1 = _iRowsIn1 * _iColsIn1;
117 int iSize2 = _iRowsIn2 * _iColsIn2;
118
119 int i = 0;
120 int j = 0;
121 int k = 0;
122 int l = 0;
123 int idx = 0;
124
125 for (i = 0 ; i < iSize1 ; i += _iRowsIn1)
126 {
127 int i2 = i + _iRowsIn1;
128 for (j = 0 ; j < iSize2 ; j += _iRowsIn2)
129 {
130 int j2 = j + _iRowsIn2;
131 for (k = i ; k < i2 ; k++)
132 {
133 //double c = _pdblDataIn1[k];
134 for (l = j ; l < j2 ; l++)
135 {
136 _pdblRealOut[idx] = _pdblRealIn1[k] * _pdblRealIn2[l] - _pdblImgIn1[k] * _pdblImgIn2[l];
137 _pdblImgOut[idx] = _pdblRealIn1[k] * _pdblImgIn2[l] + _pdblImgIn1[k] * _pdblRealIn2[l];
138 idx++;
139 }
140 }
141 }
142 }
143 }
144
145 // convert : a => 1 ./ a
conv_real_input(double * _pdblData,int _iSize)146 int conv_real_input(double* _pdblData, int _iSize)
147 {
148 int i;
149 for (i = 0 ; i < _iSize ; i++)
150 {
151 if (_pdblData[i] != 0)
152 {
153 _pdblData[i] = 1.0 / _pdblData[i];
154 }
155 else
156 {
157 return 1;
158 }
159 }
160 return 0;
161 }
162
conv_img_input(double * _pdblReal,double * _pdblImg,int _iSize)163 int conv_img_input(double* _pdblReal, double* _pdblImg, int _iSize)
164 {
165 int i;
166 for (i = 0 ; i < _iSize ; i++)
167 {
168 double dblR = _pdblReal[i];
169 double dblI = _pdblImg[i];
170
171 double dblTemp = dblR * dblR + dblI * dblI;
172 if (dblTemp != 0)
173 {
174 _pdblReal[i] = _pdblReal[i] / dblTemp;
175 _pdblImg[i] = - _pdblImg[i] / dblTemp;
176 }
177 else
178 {
179 return 1;
180 }
181 }
182 return 0;
183 }
184
185