CVM Class Library  8.1
This C++ class library encapsulates concepts of vector and different matrices including square, band, symmetric and hermitian ones in Euclidean space of real and complex numbers.
 All Classes Files Functions Variables Typedefs Friends Macros Pages
arrays.cpp
Go to the documentation of this file.
1 // CVM Class Library
2 // http://cvmlib.com
3 //
4 // Copyright Sergei Nikolaev 1992-2014
5 // Distributed under the Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
8 
9 #include "cvm.h"
10 #include "blas.h"
11 #include <cmath>
12 
13 CVM_NAMESPACE_BEG
14 
15 // 8.0: tint-based spec due to Array->basic_array refactoring
16 template<>
17 CVM_API tint __norm<tint, tint>(const tint* pd, tint nSize, tint nIncr)
18 {
19  CVM_ASSERT(pd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(tint))
20  rvector rv(nSize);
21  for (tint i = 0; i < nSize; ++i) {
22  rv[CVM0+i] = static_cast<treal>(pd[i*nIncr]);
23  }
24  return static_cast<tint>(floor(rv.norm() + 0.5));
25 }
26 
27 // 8.0: tint-based spec due to Array->basic_array refactoring
28 template<>
29 CVM_API void __scal<tint, tint>(tint* pd, tint nSize, tint nIncr, tint dScal)
30 {
31  CVM_ASSERT(pd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(tint))
32  for (tint i = 0; i < nSize; ++i) {
33  pd[i*nIncr] *= dScal;
34  }
35 }
36 
37 // 8.0: tint-based spec due to Array->basic_array refactoring
38 template<>
39 CVM_API tint __idamax<tint>(const tint* pd, tint nSize, tint nIncr)
40 {
41  CVM_ASSERT(pd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(tint))
42  tint ret = 0, max = - (std::numeric_limits<tint>::max)(), val;
43  for (tint i = 0; i < nSize; ++i) {
44  val = pd[i*nIncr];
45  if (val > max) {
46  max = val;
47  ret = i;
48  }
49  }
50  return ret + CVM0;
51 }
52 
53 // 8.0: tint-based spec due to Array->basic_array refactoring
54 template<>
55 CVM_API tint __idamin<tint>(const tint* pd, tint nSize, tint nIncr)
56 {
57  CVM_ASSERT(pd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(tint))
58  tint ret = 0, min = (std::numeric_limits<tint>::max)(), val;
59  for (tint i = 0; i < nSize; ++i) {
60  val = pd[i*nIncr];
61  if (val < min) {
62  min = val;
63  ret = i;
64  }
65  }
66  return ret + CVM0;
67 }
68 
69 
70 template<>
71 CVM_API float __norm<float, float>(const float* mpd, tint nSize, tint nIncr)
72 {
73  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(float))
74  return SNRM2 (&nSize, mpd, &nIncr);
75 }
76 
77 template<>
78 CVM_API double __norm<double, double>(const double* mpd, tint nSize, tint nIncr)
79 {
80  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(double))
81  return DNRM2 (&nSize, mpd, &nIncr);
82 }
83 
84 template<>
85 CVM_API float __norm<float, std::complex<float> >(const std::complex<float>* mpd, tint nSize, tint nIncr)
86 {
87  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<float>))
88  return SCNRM2 (&nSize, mpd, &nIncr);
89 }
90 
91 template<>
92 CVM_API double __norm<double, std::complex<double> >(const std::complex<double>* mpd, tint nSize, tint nIncr)
93 {
94  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<double>))
95  return DZNRM2 (&nSize, mpd, &nIncr);
96 }
97 
98 template<>
99 CVM_API tint __idamax<float>(const float* mpd, tint nSize, tint nIncr)
100 {
101  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(float))
102  return ISAMAX (&nSize, mpd, &nIncr) - (1 - CVM0);
103 }
104 
105 template<>
106 CVM_API tint __idamax<double>(const double* mpd, tint nSize, tint nIncr)
107 {
108  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(double))
109  return IDAMAX (&nSize, mpd, &nIncr) - (1 - CVM0);
110 }
111 
112 template<>
113 CVM_API tint __idamin<float>(const float* mpd, tint nSize, tint nIncr)
114 {
115  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(float))
116  return ISAMIN (&nSize, mpd, &nIncr) - (1 - CVM0);
117 }
118 
119 template<>
120 CVM_API tint __idamin<double>(const double* mpd, tint nSize, tint nIncr)
121 {
122  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(double))
123  return IDAMIN (&nSize, mpd, &nIncr) - (1 - CVM0);
124 }
125 
126 template<>
127 CVM_API tint __idamax<std::complex<float> >(const std::complex<float>* mpd, tint nSize, tint nIncr)
128 {
129  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<float>))
130  return ICAMAX (&nSize, mpd, &nIncr) - (1 - CVM0);
131 }
132 
133 template<>
134 CVM_API tint __idamax<std::complex<double> >(const std::complex<double>* mpd, tint nSize, tint nIncr)
135 {
136  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<double>))
137  return IZAMAX (&nSize, mpd, &nIncr) - (1 - CVM0);
138 }
139 
140 template<>
141 CVM_API tint __idamin<std::complex<float> >(const std::complex<float>* mpd, tint nSize, tint nIncr)
142 {
143  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<float>))
144  return ICAMIN (&nSize, mpd, &nIncr) - (1 - CVM0);
145 }
146 
147 template<>
148 CVM_API tint __idamin<std::complex<double> >(const std::complex<double>* mpd, tint nSize, tint nIncr)
149 {
150  CVM_ASSERT(mpd, ((nSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<double>))
151  return IZAMIN (&nSize, mpd, &nIncr) - (1 - CVM0);
152 }
153 
154 template<>
155 CVM_API void __add<float>(float* mpd, tint mnSize, tint mnIncr, const float* pv, tint nIncr)
156 {
157  static const float one(1.F);
158  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(float))
159  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(float))
160  SAXPY (&mnSize, &one, pv, &nIncr, mpd, &mnIncr);
161 }
162 
163 template<>
164 CVM_API void __add<double>(double* mpd, tint mnSize, tint mnIncr, const double* pv, tint nIncr)
165 {
166  static const double one(1.);
167  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(double))
168  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(double))
169  DAXPY (&mnSize, &one, pv, &nIncr, mpd, &mnIncr);
170 }
171 
172 template<>
173 CVM_API void __subtract<float>(float* mpd, tint mnSize, tint mnIncr, const float* pv, tint nIncr)
174 {
175  static const float mone(-1.F);
176  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(float))
177  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(float))
178  SAXPY (&mnSize, &mone, pv, &nIncr, mpd, &mnIncr);
179 }
180 
181 template<>
182 CVM_API void __subtract<double>(double* mpd, tint mnSize, tint mnIncr, const double* pv, tint nIncr)
183 {
184  static const double mone(-1.F);
185  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(double))
186  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(double))
187  DAXPY (&mnSize, &mone, pv, &nIncr, mpd, &mnIncr);
188 }
189 
190 template<>
191 CVM_API void __add<std::complex<float> >(std::complex<float>* mpd, tint mnSize, tint mnIncr, const std::complex<float>* pv, tint nIncr)
192 {
193  static const std::complex<float> one(1.F, 0.F);
194  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<float>))
195  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<float>))
196  CAXPY (&mnSize, &one, pv, &nIncr, mpd, &mnIncr);
197 }
198 
199 template<>
200 CVM_API void __add<std::complex<double> >(std::complex<double>* mpd, tint mnSize, tint mnIncr, const std::complex<double>* pv, tint nIncr)
201 {
202  static const std::complex<double> one(1., 0.);
203  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<double>))
204  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<double>))
205  ZAXPY (&mnSize, &one, pv, &nIncr, mpd, &mnIncr);
206 }
207 
208 template<>
209 CVM_API void __subtract<std::complex<float> >(std::complex<float>* mpd, tint mnSize, tint mnIncr, const std::complex<float>* pv, tint nIncr)
210 {
211  static const std::complex<float> mone(-1.F, 0.F);
212  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<float>))
213  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<float>))
214  CAXPY (&mnSize, &mone, pv, &nIncr, mpd, &mnIncr);
215 }
216 
217 template<>
218 CVM_API void __subtract<std::complex<double> >(std::complex<double>* mpd, tint mnSize, tint mnIncr, const std::complex<double>* pv, tint nIncr)
219 {
220  static const std::complex<double> mone(-1., 0.);
221  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<double>))
222  CVM_ASSERT(pv, ((mnSize - CVM0) * nIncr + CVM0) * sizeof(std::complex<double>))
223  ZAXPY (&mnSize, &mone, pv, &nIncr, mpd, &mnIncr);
224 }
225 
226 template<>
227 CVM_API void __scal<float, float>(float* mpd, tint mnSize, tint mnIncr, float dScal)
228 {
229  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(float))
230  SSCAL (&mnSize, &dScal, mpd, &mnIncr);
231 }
232 
233 template<>
234 CVM_API void __scal<double, double>(double* mpd, tint mnSize, tint mnIncr, double dScal)
235 {
236  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(double))
237  DSCAL (&mnSize, &dScal, mpd, &mnIncr);
238 }
239 
240 template<>
241 CVM_API void __scal<float, std::complex<float> >(std::complex<float>* mpd, tint mnSize, tint mnIncr, float dScal)
242 {
243  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<float>))
244  CSSCAL (&mnSize, &dScal, mpd, &mnIncr);
245 }
246 
247 template<>
248 CVM_API void __scal<double, std::complex<double> >(std::complex<double>* mpd, tint mnSize, tint mnIncr, double dScal)
249 {
250  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<double>))
251  ZDSCAL (&mnSize, &dScal, mpd, &mnIncr);
252 }
253 
255 template<>
256 CVM_API void __scal<std::complex<float>, std::complex<float> >(std::complex<float>* mpd, tint mnSize, tint mnIncr, std::complex<float> cScal)
257 {
258  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<float>))
259  CSCAL (&mnSize, &cScal, mpd, &mnIncr);
260 }
261 
262 template<>
263 CVM_API void __scal<std::complex<double>, std::complex<double> >(std::complex<double>* mpd, tint mnSize, tint mnIncr, std::complex<double> cScal)
264 {
265  CVM_ASSERT(mpd, ((mnSize - CVM0) * mnIncr + CVM0) * sizeof(std::complex<double>))
266  ZSCAL (&mnSize, &cScal, mpd, &mnIncr);
267 }
269 
270 template<>
271 CVM_API void __copy2<float, std::complex<float> > (std::complex<float>* mpd, tint mnSize, tint mnIncr,
272  const float* pRe, const float* pIm, tint nReIncr, tint nImIncr)
273 {
274  const tint nIncr2 = mnIncr * 2;
275  float* pR = __get_real_p<float>(mpd);
276  float* pI = __get_imag_p<float>(mpd);
277 
278  if (pRe)
279  {
280  CVM_ASSERT(pRe, ((mnSize - CVM0) * nReIncr + CVM0) * sizeof(float))
281  CVM_ASSERT(pR, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(float))
282  SCOPY (&mnSize, pRe, &nReIncr, pR, &nIncr2);
283  }
284  else
285  {
286  static const float zero(0.F);
287  CVM_ASSERT(pR, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(float))
288  SSCAL (&mnSize, &zero, pR, &nIncr2);
289  }
290 
291  if (pIm)
292  {
293  CVM_ASSERT(pIm, ((mnSize - CVM0) * nImIncr + CVM0) * sizeof(float))
294  CVM_ASSERT(pI, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(float))
295  SCOPY (&mnSize, pIm, &nImIncr, pI, &nIncr2);
296  }
297  else
298  {
299  static const float zero(0.F);
300  CVM_ASSERT(pI, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(float))
301  SSCAL (&mnSize, &zero, pI, &nIncr2);
302  }
303 }
304 
305 template<>
306 CVM_API void __copy2<double, std::complex<double> > (std::complex<double>* mpd, tint mnSize, tint mnIncr,
307  const double* pRe, const double* pIm, tint nReIncr, tint nImIncr)
308 {
309  const tint nIncr2 = mnIncr * 2;
310  double* pR = __get_real_p<double>(mpd);
311  double* pI = __get_imag_p<double>(mpd);
312 
313  if (pRe)
314  {
315  CVM_ASSERT(pRe, ((mnSize - CVM0) * nReIncr + CVM0) * sizeof(double))
316  CVM_ASSERT(pR, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(double))
317  DCOPY (&mnSize, pRe, &nReIncr, pR, &nIncr2);
318  }
319  else
320  {
321  static const double zero(0.);
322  CVM_ASSERT(pR, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(double))
323  DSCAL (&mnSize, &zero, pR, &nIncr2);
324  }
325 
326  if (pIm)
327  {
328  CVM_ASSERT(pIm, ((mnSize - CVM0) * nImIncr + CVM0) * sizeof(double))
329  CVM_ASSERT(pI, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(double))
330  DCOPY (&mnSize, pIm, &nImIncr, pI, &nIncr2);
331  }
332  else
333  {
334  static const double zero(0.);
335  CVM_ASSERT(pI, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(double))
336  DSCAL (&mnSize, &zero, pI, &nIncr2);
337  }
338 }
339 
340 template<>
341 CVM_API void __copy_real<float, std::complex<float> > (std::complex<float>* mpd, tint mnSize, tint mnIncr, const float* pRe, tint nReIncr)
342 {
343  float* pdr = __get_real_p<float>(mpd);
344  if (pdr != pRe)
345  {
346  const tint nIncr2 = mnIncr * 2;
347  CVM_ASSERT(pRe, ((mnSize - CVM0) * nReIncr + CVM0) * sizeof(float))
348  CVM_ASSERT(pdr, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(float))
349  SCOPY (&mnSize, pRe, &nReIncr, pdr, &nIncr2);
350  }
351 }
352 
353 template<>
354 CVM_API void __copy_real<double, std::complex<double> > (std::complex<double>* mpd, tint mnSize, tint mnIncr, const double* pRe, tint nReIncr)
355 {
356  double* pdr = __get_real_p<double>(mpd);
357  if (pdr != pRe)
358  {
359  const tint nIncr2 = mnIncr * 2;
360  CVM_ASSERT(pRe, ((mnSize - CVM0) * nReIncr + CVM0) * sizeof(double))
361  CVM_ASSERT(pdr, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(double))
362  DCOPY (&mnSize, pRe, &nReIncr, pdr, &nIncr2);
363  }
364 }
365 
366 template<>
367 CVM_API void __copy_imag<float, std::complex<float> > (std::complex<float>* mpd, tint mnSize, tint mnIncr, const float* pIm, tint nImIncr)
368 {
369  float* pdi = __get_imag_p<float>(mpd);
370  if (pdi != pIm)
371  {
372  const tint nIncr2 = mnIncr * 2;
373  CVM_ASSERT(pIm, ((mnSize - CVM0) * nImIncr + CVM0) * sizeof(float))
374  CVM_ASSERT(pdi, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(float))
375  SCOPY (&mnSize, pIm, &nImIncr, pdi, &nIncr2);
376  }
377 }
378 
379 template<>
380 CVM_API void __copy_imag<double, std::complex<double> > (std::complex<double>* mpd, tint mnSize, tint mnIncr, const double* pIm, tint nImIncr)
381 {
382  double* pdi = __get_imag_p<double>(mpd);
383  if (pdi != pIm)
384  {
385  const tint nIncr2 = mnIncr * 2;
386  CVM_ASSERT(pIm, ((mnSize - CVM0) * nImIncr + CVM0) * sizeof(double))
387  CVM_ASSERT(pdi, ((mnSize - CVM0) * nIncr2 + CVM0) * sizeof(double))
388  DCOPY (&mnSize, pIm, &nImIncr, pdi, &nIncr2);
389  }
390 }
391 
392 template<>
393 CVM_API void __conj<std::complex<float> > (std::complex<float>* mpd, tint mnSize, tint mnIncr)
394 {
395  CLACGV (&mnSize, mpd, &mnIncr);
396 }
397 
398 template<>
399 CVM_API void __conj<std::complex<double> > (std::complex<double>* mpd, tint mnSize, tint mnIncr)
400 {
401  ZLACGV (&mnSize, mpd, &mnIncr);
402 }
403 
404 template<>
405 CVM_API void __randomize<float> (float* mpd, tint mnSize, tint mnIncr, float dFrom, float dTo)
406 {
407  const tint nSize = mnSize * mnIncr;
408  for (tint i = 0; i < nSize; i += mnIncr)
409  {
410  mpd[i] = Randomizer<float>::get (dFrom, dTo);
411  }
412 }
413 
414 template<>
415 CVM_API void __randomize<double> (double* mpd, tint mnSize, tint mnIncr, double dFrom, double dTo)
416 {
417  const tint nSize = mnSize * mnIncr;
418  for (tint i = 0; i < nSize; i += mnIncr)
419  {
420  mpd[i] = Randomizer<double>::get (dFrom, dTo);
421  }
422 }
423 
424 template<>
425 CVM_API void __randomize_real<std::complex<float>, float> (std::complex<float>* mpd, tint mnSize, tint mnIncr,
426  float dFrom, float dTo)
427 {
428  const tint nIncr = 2 * mnIncr;
429  const tint nSize = mnSize * nIncr;
430  float* pdr = __get_real_p<float>(mpd);
431 
432  for (tint i = 0; i < nSize; i += nIncr)
433  {
434  pdr[i] = Randomizer<float>::get (dFrom, dTo);
435  }
436 }
437 
438 template<>
439 CVM_API void __randomize_real<std::complex<double>, double> (std::complex<double>* mpd, tint mnSize, tint mnIncr,
440  double dFrom, double dTo)
441 {
442  const tint nIncr = 2 * mnIncr;
443  const tint nSize = mnSize * nIncr;
444  double* pdr = __get_real_p<double>(mpd);
445 
446  for (tint i = 0; i < nSize; i += nIncr)
447  {
448  pdr[i] = Randomizer<double>::get (dFrom, dTo);
449  }
450 }
451 
452 template<>
453 CVM_API void __randomize_imag<std::complex<float>, float> (std::complex<float>* mpd, tint mnSize, tint mnIncr,
454  float dFrom, float dTo)
455 {
456  const tint nIncr = 2 * mnIncr;
457  const tint nSize = mnSize * nIncr;
458  float* pdi = __get_imag_p<float>(mpd);
459 
460  for (tint i = 0; i < nSize; i += nIncr)
461  {
462  pdi[i] = Randomizer<float>::get (dFrom, dTo);
463  }
464 }
465 
466 template<>
467 CVM_API void __randomize_imag<std::complex<double>, double> (std::complex<double>* mpd, tint mnSize, tint mnIncr,
468  double dFrom, double dTo)
469 {
470  const tint nIncr = 2 * mnIncr;
471  const tint nSize = mnSize * nIncr;
472  double* pdi = __get_imag_p<double>(mpd);
473 
474  for (tint i = 0; i < nSize; i += nIncr)
475  {
476  pdi[i] = Randomizer<double>::get (dFrom, dTo);
477  }
478 }
479 
480 CVM_NAMESPACE_END
481 
Generated on Tue Jul 22 2014 18:52:07 for CVM Class Library by   doxygen 1.8.1.2