73 complex(
const double &_x,
const double &_y):
x(_x),
y(_y){};
91 if( fabs(z.y)<fabs(z.x) )
112const complex
operator/(
const complex& lhs,
const complex& rhs);
113bool operator==(
const complex& lhs,
const complex& rhs);
114bool operator!=(
const complex& lhs,
const complex& rhs);
115const complex
operator+(
const complex& lhs);
116const complex
operator-(
const complex& lhs);
117const complex
operator+(
const complex& lhs,
const complex& rhs);
118const complex
operator+(
const complex& lhs,
const double& rhs);
119const complex
operator+(
const double& lhs,
const complex& rhs);
120const complex
operator-(
const complex& lhs,
const complex& rhs);
121const complex
operator-(
const complex& lhs,
const double& rhs);
122const complex
operator-(
const double& lhs,
const complex& rhs);
123const complex
operator*(
const complex& lhs,
const complex& rhs);
124const complex
operator*(
const complex& lhs,
const double& rhs);
125const complex
operator*(
const double& lhs,
const complex& rhs);
126const complex
operator/(
const complex& lhs,
const complex& rhs);
127const complex
operator/(
const double& lhs,
const complex& rhs);
128const complex
operator/(
const complex& lhs,
const double& rhs);
130const complex
conj(
const complex &z);
131const complex
csqr(
const complex &z);
145class const_raw_vector
176class raw_vector :
public const_raw_vector<T>
190T vdotproduct(const_raw_vector<T> v1, const_raw_vector<T> v2)
193 if( v1.GetStep()==1 && v2.GetStep()==1 )
199 const T *p1 = v1.GetData();
200 const T *p2 = v2.GetData();
201 int imax = v1.GetLength()/4;
203 for(
i=imax;
i!=0;
i--)
205 r += (*p1)*(*p2) + p1[1]*p2[1] + p1[2]*p2[2] + p1[3]*p2[3];
209 for(
i=0;
i<v1.GetLength()%4;
i++)
210 r += (*(p1++))*(*(p2++));
218 int offset11 = v1.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
219 int offset21 = v2.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
221 const T *p1 = v1.GetData();
222 const T *p2 = v2.GetData();
223 int imax = v1.GetLength()/4;
225 for(
i=0;
i<imax;
i++)
227 r += (*p1)*(*p2) + p1[offset11]*p2[offset21] + p1[offset12]*p2[offset22] + p1[offset13]*p2[offset23];
231 for(
i=0;
i<v1.GetLength()%4;
i++)
249 if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
254 T *p1 = vdst.GetData();
255 const T *p2 = vsrc.GetData();
256 int imax = vdst.GetLength()/2;
258 for(
i=imax;
i!=0;
i--)
265 if(vdst.GetLength()%2 != 0)
274 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
275 int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
276 T *p1 = vdst.GetData();
277 const T *p2 = vsrc.GetData();
278 int imax = vdst.GetLength()/4;
280 for(
i=0;
i<imax;
i++)
283 p1[offset11] = p2[offset21];
284 p1[offset12] = p2[offset22];
285 p1[offset13] = p2[offset23];
289 for(
i=0;
i<vdst.GetLength()%4;
i++)
292 p1 += vdst.GetStep();
293 p2 += vsrc.GetStep();
307 if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
312 T *p1 = vdst.GetData();
313 const T *p2 = vsrc.GetData();
314 int imax = vdst.GetLength()/2;
316 for(
i=0;
i<imax;
i++)
323 if(vdst.GetLength()%2 != 0)
332 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
333 int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
334 T *p1 = vdst.GetData();
335 const T *p2 = vsrc.GetData();
336 int imax = vdst.GetLength()/4;
338 for(
i=imax;
i!=0;
i--)
341 p1[offset11] = -p2[offset21];
342 p1[offset12] = -p2[offset22];
343 p1[offset13] = -p2[offset23];
347 for(
i=0;
i<vdst.GetLength()%4;
i++)
350 p1 += vdst.GetStep();
351 p2 += vsrc.GetStep();
361template<
class T,
class T2>
365 if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
370 T *p1 = vdst.GetData();
371 const T *p2 = vsrc.GetData();
372 int imax = vdst.GetLength()/4;
374 for(
i=imax;
i!=0;
i--)
383 for(
i=0;
i<vdst.GetLength()%4;
i++)
384 *(p1++) =
alpha*(*(p2++));
392 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
393 int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
394 T *p1 = vdst.GetData();
395 const T *p2 = vsrc.GetData();
396 int imax = vdst.GetLength()/4;
398 for(
i=0;
i<imax;
i++)
401 p1[offset11] =
alpha*p2[offset21];
402 p1[offset12] =
alpha*p2[offset22];
403 p1[offset13] =
alpha*p2[offset23];
407 for(
i=0;
i<vdst.GetLength()%4;
i++)
410 p1 += vdst.GetStep();
411 p2 += vsrc.GetStep();
425 if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
430 T *p1 = vdst.GetData();
431 const T *p2 = vsrc.GetData();
432 int imax = vdst.GetLength()/4;
434 for(
i=imax;
i!=0;
i--)
443 for(
i=0;
i<vdst.GetLength()%4;
i++)
452 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
453 int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
454 T *p1 = vdst.GetData();
455 const T *p2 = vsrc.GetData();
456 int imax = vdst.GetLength()/4;
458 for(
i=0;
i<imax;
i++)
461 p1[offset11] += p2[offset21];
462 p1[offset12] += p2[offset22];
463 p1[offset13] += p2[offset23];
467 for(
i=0;
i<vdst.GetLength()%4;
i++)
470 p1 += vdst.GetStep();
471 p2 += vsrc.GetStep();
481template<
class T,
class T2>
485 if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
490 T *p1 = vdst.GetData();
491 const T *p2 = vsrc.GetData();
492 int imax = vdst.GetLength()/4;
494 for(
i=imax;
i!=0;
i--)
497 p1[1] +=
alpha*p2[1];
498 p1[2] +=
alpha*p2[2];
499 p1[3] +=
alpha*p2[3];
503 for(
i=0;
i<vdst.GetLength()%4;
i++)
504 *(p1++) +=
alpha*(*(p2++));
512 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
513 int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
514 T *p1 = vdst.GetData();
515 const T *p2 = vsrc.GetData();
516 int imax = vdst.GetLength()/4;
518 for(
i=0;
i<imax;
i++)
521 p1[offset11] +=
alpha*p2[offset21];
522 p1[offset12] +=
alpha*p2[offset22];
523 p1[offset13] +=
alpha*p2[offset23];
527 for(
i=0;
i<vdst.GetLength()%4;
i++)
530 p1 += vdst.GetStep();
531 p2 += vsrc.GetStep();
545 if( vdst.GetStep()==1 && vsrc.GetStep()==1 )
550 T *p1 = vdst.GetData();
551 const T *p2 = vsrc.GetData();
552 int imax = vdst.GetLength()/4;
554 for(
i=imax;
i!=0;
i--)
563 for(
i=0;
i<vdst.GetLength()%4;
i++)
572 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
573 int offset21 = vsrc.GetStep(), offset22 = 2*offset21, offset23 = 3*offset21, offset24 = 4*offset21;
574 T *p1 = vdst.GetData();
575 const T *p2 = vsrc.GetData();
576 int imax = vdst.GetLength()/4;
578 for(
i=0;
i<imax;
i++)
581 p1[offset11] -= p2[offset21];
582 p1[offset12] -= p2[offset22];
583 p1[offset13] -= p2[offset23];
587 for(
i=0;
i<vdst.GetLength()%4;
i++)
590 p1 += vdst.GetStep();
591 p2 += vsrc.GetStep();
601template<
class T,
class T2>
611template<
class T,
class T2>
614 if( vdst.GetStep()==1 )
619 T *p1 = vdst.GetData();
620 int imax = vdst.GetLength()/4;
622 for(
i=imax;
i!=0;
i--)
630 for(
i=0;
i<vdst.GetLength()%4;
i++)
639 int offset11 = vdst.GetStep(), offset12 = 2*offset11, offset13 = 3*offset11, offset14 = 4*offset11;
640 T *p1 = vdst.GetData();
641 int imax = vdst.GetLength()/4;
643 for(
i=0;
i<imax;
i++)
646 p1[offset11] *=
alpha;
647 p1[offset12] *=
alpha;
648 p1[offset13] *=
alpha;
651 for(
i=0;
i<vdst.GetLength()%4;
i++)
654 p1 += vdst.GetStep();
688 #ifndef UNSAFE_MEM_COPY
713 #ifndef UNSAFE_MEM_COPY
758 for(
int i=iLow;
i<=iHigh;
i++)
759 (*
this)(
i) = pContent[
i-iLow];
815class template_2d_array
842 #ifndef UNSAFE_MEM_COPY
869 #ifndef UNSAFE_MEM_COPY
899 void setbounds(
int iLow1,
int iHigh1,
int iLow2,
int iHigh2 )
903 m_iVecSize = (iHigh1-iLow1+1)*(iHigh2-iLow2+1);
913 void setcontent(
int iLow1,
int iHigh1,
int iLow2,
int iHigh2,
const T *pContent )
953 return raw_vector<T>(&((*
this)(iRow, iColumnStart)), iColumnEnd-iColumnStart+1, 1);
1007double sqr(
double x);
1008int maxint(
int m1,
int m2);
1009int minint(
int m1,
int m2);
1010double maxreal(
double m1,
double m2);
1011double minreal(
double m1,
double m2);
1095 template<
unsigned int Precision>
1105 if(
rval->refCount==0 )
1142#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
1143 template<
unsigned int Precision2>
1176 if(
rval->refCount==0 )
1183#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
1184 template<
unsigned int Precision2>
1187 if( (
void*)
this==(
void*)(&r) )
1258 template<
unsigned int Precision>
1263 WerrorS(
"incorrectPrecision");
1266 template<
unsigned int Precision>
1271 mpfr_set_ui(getWritePtr(), 0, GMP_RNDN);
1274 template<
unsigned int Precision>
1279 mpfr_set_si(getWritePtr(), sv, GMP_RNDN);
1282 template<
unsigned int Precision>
1287 mpfr_set_ui(getWritePtr(),
v, GMP_RNDN);
1290 template<
unsigned int Precision>
1295 mpfr_set_ld(getWritePtr(),
v, GMP_RNDN);
1298 template<
unsigned int Precision>
1303 mpfr_strtofr(getWritePtr(),
s,
NULL, 0, GMP_RNDN);
1306 template<
unsigned int Precision>
1312 template<
unsigned int Precision>
1315 if( rval->refCount==1 )
1318 mpfr_set(newrval->value, rval->value, GMP_RNDN);
1324 template<
unsigned int Precision>
1327 return mpfr_number_p(getReadPtr())!=0;
1330 template<
unsigned int Precision>
1333 if( !isFiniteNumber() )
1335 return mpfr_sgn(getReadPtr())>0;
1338 template<
unsigned int Precision>
1341 return mpfr_zero_p(getReadPtr())!=0;
1344 template<
unsigned int Precision>
1347 if( !isFiniteNumber() )
1349 return mpfr_sgn(getReadPtr())<0;
1352 template<
unsigned int Precision>
1355 return getUlpOf(*
this);
1358 template<
unsigned int Precision>
1361 return mpfr_get_d(getReadPtr(), GMP_RNDN);
1364 template<
unsigned int Precision>
1370 if( !isFiniteNumber() )
1375 ptr = mpfr_get_str(
NULL, &_e, 16, 0, getReadPtr(), GMP_RNDN);
1386 signed long iexpval;
1390 ptr = mpfr_get_str(
NULL, &expval, 16, 0, getReadPtr(), GMP_RNDN);
1393 if( iexpval!=expval )
1396 sprintf(buf_e,
"%ld",
long(iexpval));
1406 mpfr_free_str(ptr2);
1410 template<
unsigned int Precision>
1418 if( !isFiniteNumber() )
1423 ptr = mpfr_get_str(
NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
1434 signed long iexpval;
1438 ptr = mpfr_get_str(
NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
1441 if( iexpval!=expval )
1444 sprintf(buf_e,
"%ld",
long(iexpval));
1454 mpfr_free_str(ptr2);
1457 template<
unsigned int Precision>
1460 char *toString_Block=(
char *)
omAlloc(256);
1464 if( !isFiniteNumber() )
1468 ptr = mpfr_get_str(
NULL, &_e, 10, 0, getReadPtr(), GMP_RNDN);
1469 strcpy(toString_Block, ptr);
1471 return toString_Block;
1479 signed long iexpval;
1483 ptr = mpfr_get_str(
NULL, &expval, 10, 0, getReadPtr(), GMP_RNDN);
1486 if( iexpval!=expval )
1489 sprintf(buf_e,
"%ld",
long(iexpval));
1493 sprintf(toString_Block,
"-0.%sE%s",ptr,buf_e);
1496 sprintf(toString_Block,
"0.%sE%s",ptr,buf_e);
1497 mpfr_free_str(ptr2);
1498 return toString_Block;
1501 template<
unsigned int Precision>
1504 if( !
x.isFiniteNumber() )
1509 mpfr_nextabove(r.getWritePtr());
1510 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1514 mpfr_get_exp(
x.getReadPtr()),
1524 template<
unsigned int Precision>
1528 mpfr_nextabove(r.getWritePtr());
1529 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1533 template<
unsigned int Precision>
1537 mpfr_nextabove(r.getWritePtr());
1538 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1547 template<
unsigned int Precision>
1551 mpfr_nextabove(r.getWritePtr());
1552 mpfr_sub_ui(r.getWritePtr(), r.getWritePtr(), 1, GMP_RNDN);
1561 template<
unsigned int Precision>
1565 mpfr_nextbelow(r.getWritePtr());
1566 mpfr_set_exp(r.getWritePtr(),mpfr_get_emax());
1570 template<
unsigned int Precision>
1574 mpfr_set_exp(r.getWritePtr(),mpfr_get_emin());
1578 template<
unsigned int Precision>
1584 template<
unsigned int Precision>
1588 mp_exp_t e1 = mpfr_get_emax();
1589 mp_exp_t e2 = -mpfr_get_emin();
1590 mp_exp_t e = e1>e2 ? e1 : e2;
1591 mpfr_set_exp(r.getWritePtr(), e-5);
1595 template<
unsigned int Precision>
1599 mp_exp_t e1 = mpfr_get_emax();
1600 mp_exp_t e2 = -mpfr_get_emin();
1601 mp_exp_t e = e1>e2 ? e1 : e2;
1602 mpfr_set_exp(r.getWritePtr(), 2-(e-5));
1606 template<
unsigned int Precision>
1617 template<
unsigned int Precision>
1623 template<
unsigned int Precision>
1629 template<
unsigned int Precision>
1635 template<
unsigned int Precision>
1641 template<
unsigned int Precision>
1647 template<
unsigned int Precision>
1656 template<
unsigned int Precision>
1657 const ampf<Precision>
operator+(
const ampf<Precision>& op1)
1662 template<
unsigned int Precision>
1666 mpfr_neg(
v->value, op1.getReadPtr(), GMP_RNDN);
1670 template<
unsigned int Precision>
1674 mpfr_add(
v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1678 template<
unsigned int Precision>
1682 mpfr_sub(
v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1687 template<
unsigned int Precision>
1691 mpfr_mul(
v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1695 template<
unsigned int Precision>
1699 mpfr_div(
v->value, op1.getReadPtr(), op2.getReadPtr(), GMP_RNDN);
1706 template<
unsigned int Precision>
1711 mpfr_sqr(
res.getWritePtr(),
x.getReadPtr(), GMP_RNDN);
1715 template<
unsigned int Precision>
1718 int s = mpfr_sgn(
x.getReadPtr());
1726 template<
unsigned int Precision>
1727 const ampf<Precision>
abs(
const ampf<Precision> &
x)
1730 ampf<Precision>
res;
1731 mpfr_abs(
res.getWritePtr(),
x.getReadPtr(), GMP_RNDN);
1735 template<
unsigned int Precision>
1740 mpfr_max(
res.getWritePtr(),
x.getReadPtr(),
y.getReadPtr(), GMP_RNDN);
1744 template<
unsigned int Precision>
1749 mpfr_min(
res.getWritePtr(),
x.getReadPtr(),
y.getReadPtr(), GMP_RNDN);
1753 template<
unsigned int Precision>
1758 mpfr_sqrt(
res.getWritePtr(),
x.getReadPtr(), GMP_RNDN);
1762 template<
unsigned int Precision>
1771 mpfr_clear_erangeflag();
1773 if( mpfr_erangeflag_p()!=0 )
1779 template<
unsigned int Precision>
1780 const ampf<Precision>
frac(
const ampf<Precision> &
x)
1784 mpfr_frac(r.getWritePtr(),
x.getReadPtr(), GMP_RNDN);
1788 template<
unsigned int Precision>
1797 mpfr_clear_erangeflag();
1799 if( mpfr_erangeflag_p()!=0 )
1805 template<
unsigned int Precision>
1814 mpfr_clear_erangeflag();
1816 if( mpfr_erangeflag_p()!=0 )
1822 template<
unsigned int Precision>
1831 mpfr_clear_erangeflag();
1833 if( mpfr_erangeflag_p()!=0 )
1839 template<
unsigned int Precision>
1840 const ampf<Precision>
frexp2(
const ampf<Precision> &
x, mp_exp_t *
exponent)
1844 if( !
x.isFiniteNumber() )
1854 *
exponent = mpfr_get_exp(r.getReadPtr());
1855 mpfr_set_exp(r.getWritePtr(),0);
1859 template<
unsigned int Precision>
1864 mpfr_mul_2si(r.getWritePtr(),
x.getReadPtr(),
exponent, GMP_RNDN);
1871 #define __AMP_BINARY_OPI(type) \
1872 template<unsigned int Precision> const ampf<Precision> operator+(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
1873 template<unsigned int Precision> const ampf<Precision> operator+(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
1874 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const signed type& op2) { return op1+ampf<Precision>(op2); } \
1875 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const unsigned type& op2) { return op1+ampf<Precision>(op2); } \
1876 template<unsigned int Precision> const ampf<Precision> operator-(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
1877 template<unsigned int Precision> const ampf<Precision> operator-(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
1878 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const signed type& op2) { return op1-ampf<Precision>(op2); } \
1879 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const unsigned type& op2) { return op1-ampf<Precision>(op2); } \
1880 template<unsigned int Precision> const ampf<Precision> operator*(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
1881 template<unsigned int Precision> const ampf<Precision> operator*(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
1882 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const signed type& op2) { return op1*ampf<Precision>(op2); } \
1883 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const unsigned type& op2) { return op1*ampf<Precision>(op2); } \
1884 template<unsigned int Precision> const ampf<Precision> operator/(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
1885 template<unsigned int Precision> const ampf<Precision> operator/(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
1886 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const signed type& op2) { return op1/ampf<Precision>(op2); } \
1887 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const unsigned type& op2) { return op1/ampf<Precision>(op2); } \
1888 template<unsigned int Precision> bool operator==(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
1889 template<unsigned int Precision> bool operator==(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
1890 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const signed type& op2) { return op1==ampf<Precision>(op2); } \
1891 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const unsigned type& op2) { return op1==ampf<Precision>(op2); } \
1892 template<unsigned int Precision> bool operator!=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
1893 template<unsigned int Precision> bool operator!=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
1894 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const signed type& op2) { return op1!=ampf<Precision>(op2); } \
1895 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const unsigned type& op2) { return op1!=ampf<Precision>(op2); } \
1896 template<unsigned int Precision> bool operator<=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
1897 template<unsigned int Precision> bool operator<=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
1898 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const signed type& op2) { return op1<=ampf<Precision>(op2); } \
1899 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const unsigned type& op2) { return op1<=ampf<Precision>(op2); } \
1900 template<unsigned int Precision> bool operator>=(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
1901 template<unsigned int Precision> bool operator>=(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
1902 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const signed type& op2) { return op1>=ampf<Precision>(op2); } \
1903 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const unsigned type& op2) { return op1>=ampf<Precision>(op2); } \
1904 template<unsigned int Precision> bool operator<(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
1905 template<unsigned int Precision> bool operator<(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
1906 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const signed type& op2) { return op1<ampf<Precision>(op2); } \
1907 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const unsigned type& op2) { return op1<ampf<Precision>(op2); } \
1908 template<unsigned int Precision> bool operator>(const signed type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
1909 template<unsigned int Precision> bool operator>(const unsigned type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
1910 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const signed type& op2) { return op1>ampf<Precision>(op2); } \
1911 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const unsigned type& op2) { return op1>ampf<Precision>(op2); }
1916 #undef __AMP_BINARY_OPI
1917 #define __AMP_BINARY_OPF(type) \
1918 template<unsigned int Precision> const ampf<Precision> operator+(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)+op2; } \
1919 template<unsigned int Precision> const ampf<Precision> operator+(const ampf<Precision>& op1, const type& op2) { return op1+ampf<Precision>(op2); } \
1920 template<unsigned int Precision> const ampf<Precision> operator-(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)-op2; } \
1921 template<unsigned int Precision> const ampf<Precision> operator-(const ampf<Precision>& op1, const type& op2) { return op1-ampf<Precision>(op2); } \
1922 template<unsigned int Precision> const ampf<Precision> operator*(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)*op2; } \
1923 template<unsigned int Precision> const ampf<Precision> operator*(const ampf<Precision>& op1, const type& op2) { return op1*ampf<Precision>(op2); } \
1924 template<unsigned int Precision> const ampf<Precision> operator/(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)/op2; } \
1925 template<unsigned int Precision> const ampf<Precision> operator/(const ampf<Precision>& op1, const type& op2) { return op1/ampf<Precision>(op2); } \
1926 template<unsigned int Precision> bool operator==(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)==op2; } \
1927 template<unsigned int Precision> bool operator==(const ampf<Precision>& op1, const type& op2) { return op1==ampf<Precision>(op2); } \
1928 template<unsigned int Precision> bool operator!=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)!=op2; } \
1929 template<unsigned int Precision> bool operator!=(const ampf<Precision>& op1, const type& op2) { return op1!=ampf<Precision>(op2); } \
1930 template<unsigned int Precision> bool operator<=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<=op2; } \
1931 template<unsigned int Precision> bool operator<=(const ampf<Precision>& op1, const type& op2) { return op1<=ampf<Precision>(op2); } \
1932 template<unsigned int Precision> bool operator>=(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>=op2; } \
1933 template<unsigned int Precision> bool operator>=(const ampf<Precision>& op1, const type& op2) { return op1>=ampf<Precision>(op2); } \
1934 template<unsigned int Precision> bool operator<(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)<op2; } \
1935 template<unsigned int Precision> bool operator<(const ampf<Precision>& op1, const type& op2) { return op1<ampf<Precision>(op2); } \
1936 template<unsigned int Precision> bool operator>(const type& op1, const ampf<Precision>& op2) { return ampf<Precision>(op1)>op2; } \
1937 template<unsigned int Precision> bool operator>(const ampf<Precision>& op1, const type& op2) { return op1>ampf<Precision>(op2); }
1941 #undef __AMP_BINARY_OPF
1946 template<
unsigned int Precision>
1947 const ampf<Precision>
pi()
1950 mpfr_const_pi(
v->value, GMP_RNDN);
1954 template<
unsigned int Precision>
1958 mpfr_const_pi(
v->value, GMP_RNDN);
1959 mpfr_mul_2si(
v->value,
v->value, -1, GMP_RNDN);
1963 template<
unsigned int Precision>
1967 mpfr_const_pi(
v->value, GMP_RNDN);
1968 mpfr_mul_2si(
v->value,
v->value, +1, GMP_RNDN);
1972 template<
unsigned int Precision>
1976 mpfr_sin(
v->value,
x.getReadPtr(), GMP_RNDN);
1980 template<
unsigned int Precision>
1984 mpfr_cos(
v->value,
x.getReadPtr(), GMP_RNDN);
1988 template<
unsigned int Precision>
1992 mpfr_tan(
v->value,
x.getReadPtr(), GMP_RNDN);
1996 template<
unsigned int Precision>
2000 mpfr_asin(
v->value,
x.getReadPtr(), GMP_RNDN);
2004 template<
unsigned int Precision>
2008 mpfr_acos(
v->value,
x.getReadPtr(), GMP_RNDN);
2012 template<
unsigned int Precision>
2016 mpfr_atan(
v->value,
x.getReadPtr(), GMP_RNDN);
2020 template<
unsigned int Precision>
2024 mpfr_atan2(
v->value,
y.getReadPtr(),
x.getReadPtr(), GMP_RNDN);
2028 template<
unsigned int Precision>
2032 mpfr_log(
v->value,
x.getReadPtr(), GMP_RNDN);
2036 template<
unsigned int Precision>
2040 mpfr_log2(
v->value,
x.getReadPtr(), GMP_RNDN);
2044 template<
unsigned int Precision>
2048 mpfr_log10(
v->value,
x.getReadPtr(), GMP_RNDN);
2052 template<
unsigned int Precision>
2056 mpfr_exp(
v->value,
x.getReadPtr(), GMP_RNDN);
2060 template<
unsigned int Precision>
2064 mpfr_sinh(
v->value,
x.getReadPtr(), GMP_RNDN);
2068 template<
unsigned int Precision>
2072 mpfr_cosh(
v->value,
x.getReadPtr(), GMP_RNDN);
2076 template<
unsigned int Precision>
2080 mpfr_tanh(
v->value,
x.getReadPtr(), GMP_RNDN);
2084 template<
unsigned int Precision>
2088 mpfr_pow(
v->value,
x.getReadPtr(),
y.getReadPtr(), GMP_RNDN);
2095 template<
unsigned int Precision>
2114#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
2115 template<
unsigned int Prec2>
2138#ifndef _AMP_NO_TEMPLATE_CONSTRUCTORS
2139 template<
unsigned int Precision2>
2154 template<
unsigned int Precision>
2156 {
return lhs.x==rhs.x && lhs.y==rhs.y; }
2158 template<
unsigned int Precision>
2160 {
return lhs.x!=rhs.x || lhs.y!=rhs.y; }
2162 template<
unsigned int Precision>
2163 const campf<Precision>
operator+(
const campf<Precision>& lhs)
2166 template<
unsigned int Precision>
2168 { lhs.x += rhs.x; lhs.y += rhs.y;
return lhs; }
2170 template<
unsigned int Precision>
2174 template<
unsigned int Precision>
2178 template<
unsigned int Precision>
2180 { lhs.x -= rhs.x; lhs.y -= rhs.y;
return lhs; }
2182 template<
unsigned int Precision>
2186 template<
unsigned int Precision>
2189 ampf<Precision> xx(lhs.x*rhs.x), yy(lhs.y*rhs.y), mm((lhs.x+lhs.y)*(rhs.x+rhs.y));
2195 template<
unsigned int Precision>
2199 template<
unsigned int Precision>
2205 if(
abs(rhs.y)<
abs(rhs.x) )
2217 result.y = (-lhs.x+lhs.y*e)/
f;
2222 template<
unsigned int Precision>
2229 template<
unsigned int Precision>
2236 w = xabs>yabs ? xabs : yabs;
2237 v = xabs<yabs ? xabs : yabs;
2247 template<
unsigned int Precision>
2253 template<
unsigned int Precision>
2263 #define __AMP_BINARY_OPI(type) \
2264 template<unsigned int Precision> const campf<Precision> operator+ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
2265 template<unsigned int Precision> const campf<Precision> operator+ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
2266 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
2267 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
2268 template<unsigned int Precision> const campf<Precision> operator- (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
2269 template<unsigned int Precision> const campf<Precision> operator- (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
2270 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
2271 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
2272 template<unsigned int Precision> const campf<Precision> operator* (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
2273 template<unsigned int Precision> const campf<Precision> operator* (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
2274 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
2275 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
2276 template<unsigned int Precision> const campf<Precision> operator/ (const signed type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \
2277 template<unsigned int Precision> const campf<Precision> operator/ (const unsigned type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \
2278 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const signed type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
2279 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const unsigned type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
2280 template<unsigned int Precision> bool operator==(const signed type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
2281 template<unsigned int Precision> bool operator==(const unsigned type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
2282 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const signed type& op2) { return op1.x==op2 && op1.y==0; } \
2283 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const unsigned type& op2) { return op1.x==op2 && op1.y==0; } \
2284 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const signed type& op2) { return op1.x!=op2 || op1.y!=0; } \
2285 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const unsigned type& op2) { return op1.x!=op2 || op1.y!=0; } \
2286 template<unsigned int Precision> bool operator!=(const signed type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \
2287 template<unsigned int Precision> bool operator!=(const unsigned type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; }
2292 #undef __AMP_BINARY_OPI
2293 #define __AMP_BINARY_OPF(type) \
2294 template<unsigned int Precision> const campf<Precision> operator+ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1+op2.x, op2.y); } \
2295 template<unsigned int Precision> const campf<Precision> operator+ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x+op2, op1.y); } \
2296 template<unsigned int Precision> const campf<Precision> operator- (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1-op2.x, -op2.y); } \
2297 template<unsigned int Precision> const campf<Precision> operator- (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x-op2, op1.y); } \
2298 template<unsigned int Precision> const campf<Precision> operator* (const type& op1, const campf<Precision>& op2) { return campf<Precision>(op1*op2.x, op1*op2.y); } \
2299 template<unsigned int Precision> const campf<Precision> operator* (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op2*op1.x, op2*op1.y); } \
2300 template<unsigned int Precision> const campf<Precision> operator/ (const type& op1, const campf<Precision>& op2) { return campf<Precision>(ampf<Precision>(op1),ampf<Precision>(0))/op2; } \
2301 template<unsigned int Precision> const campf<Precision> operator/ (const campf<Precision>& op1, const type& op2) { return campf<Precision>(op1.x/op2, op1.y/op2); } \
2302 template<unsigned int Precision> bool operator==(const type& op1, const campf<Precision>& op2) { return op1==op2.x && op2.y==0; } \
2303 template<unsigned int Precision> bool operator==(const campf<Precision>& op1, const type& op2) { return op1.x==op2 && op1.y==0; } \
2304 template<unsigned int Precision> bool operator!=(const type& op1, const campf<Precision>& op2) { return op1!=op2.x || op2.y!=0; } \
2305 template<unsigned int Precision> bool operator!=(const campf<Precision>& op1, const type& op2) { return op1.x!=op2 || op1.y!=0; }
2310 #undef __AMP_BINARY_OPF
2315 template<
unsigned int Precision>
2319 int i, cnt = v1.GetLength();
2328 mpfr_set_ui(r->value, 0, GMP_RNDN);
2329 for(
i=0;
i<cnt;
i++)
2331 mpfr_mul(t->value, p1->getReadPtr(), p2->getReadPtr(), GMP_RNDN);
2332 mpfr_add(r->value, r->value, t->value, GMP_RNDN);
2349 template<
unsigned int Precision>
2353 int i, cnt = vDst.GetLength();
2358 for(
i=0;
i<cnt;
i++)
2361 pDst += vDst.GetStep();
2362 pSrc += vSrc.GetStep();
2366 template<
unsigned int Precision>
2370 int i, cnt = vDst.GetLength();
2373 for(
i=0;
i<cnt;
i++)
2376 mpfr_ptr
v = pDst->getWritePtr();
2377 mpfr_neg(
v,
v, GMP_RNDN);
2378 pDst += vDst.GetStep();
2379 pSrc += vSrc.GetStep();
2383 template<
unsigned int Precision,
class T2>
2387 int i, cnt = vDst.GetLength();
2391 for(
i=0;
i<cnt;
i++)
2394 mpfr_ptr
v = pDst->getWritePtr();
2395 mpfr_mul(
v,
v, a.getReadPtr(), GMP_RNDN);
2396 pDst += vDst.GetStep();
2397 pSrc += vSrc.GetStep();
2401 template<
unsigned int Precision>
2405 int i, cnt = vDst.GetLength();
2408 for(
i=0;
i<cnt;
i++)
2410 mpfr_ptr
v = pDst->getWritePtr();
2411 mpfr_srcptr vs = pSrc->getReadPtr();
2412 mpfr_add(
v,
v, vs, GMP_RNDN);
2413 pDst += vDst.GetStep();
2414 pSrc += vSrc.GetStep();
2418 template<
unsigned int Precision,
class T2>
2422 int i, cnt = vDst.GetLength();
2426 for(
i=0;
i<cnt;
i++)
2428 mpfr_ptr
v = pDst->getWritePtr();
2429 mpfr_srcptr vs = pSrc->getReadPtr();
2430 mpfr_mul(tmp.getWritePtr(), a.getReadPtr(), vs, GMP_RNDN);
2431 mpfr_add(
v,
v, tmp.getWritePtr(), GMP_RNDN);
2432 pDst += vDst.GetStep();
2433 pSrc += vSrc.GetStep();
2437 template<
unsigned int Precision>
2441 int i, cnt = vDst.GetLength();
2444 for(
i=0;
i<cnt;
i++)
2446 mpfr_ptr
v = pDst->getWritePtr();
2447 mpfr_srcptr vs = pSrc->getReadPtr();
2448 mpfr_sub(
v,
v, vs, GMP_RNDN);
2449 pDst += vDst.GetStep();
2450 pSrc += vSrc.GetStep();
2454 template<
unsigned int Precision,
class T2>
2460 template<
unsigned int Precision,
class T2>
2463 int i, cnt = vDst.GetLength();
2466 for(
i=0;
i<cnt;
i++)
2468 mpfr_ptr
v = pDst->getWritePtr();
2469 mpfr_mul(
v, a.getReadPtr(),
v, GMP_RNDN);
2470 pDst += vDst.GetStep();
2517 template<
unsigned int Precision>
2520 amp::ampf<Precision>&
tau);
2521 template<
unsigned int Precision>
2523 amp::ampf<Precision>
tau,
2524 const ap::template_1d_array< amp::ampf<Precision> >&
v,
2529 ap::template_1d_array< amp::ampf<Precision> >& work);
2530 template<
unsigned int Precision>
2532 amp::ampf<Precision>
tau,
2533 const ap::template_1d_array< amp::ampf<Precision> >&
v,
2538 ap::template_1d_array< amp::ampf<Precision> >& work);
2581 template<
unsigned int Precision>
2584 amp::ampf<Precision>&
tau)
2587 amp::ampf<Precision>
alpha;
2588 amp::ampf<Precision> xnorm;
2589 amp::ampf<Precision>
v;
2590 amp::ampf<Precision>
beta;
2591 amp::ampf<Precision> mx;
2676 template<
unsigned int Precision>
2678 amp::ampf<Precision>
tau,
2679 const ap::template_1d_array< amp::ampf<Precision> >&
v,
2684 ap::template_1d_array< amp::ampf<Precision> >& work)
2686 amp::ampf<Precision> t;
2689 if(
tau==0 || n1>n2 || m1>m2 )
2697 for(
i=n1;
i<=n2;
i++)
2701 for(
i=m1;
i<=m2;
i++)
2704 ap::vadd(work.getvector(n1, n2), c.getrow(
i, n1, n2), t);
2710 for(
i=m1;
i<=m2;
i++)
2713 ap::vsub(c.getrow(
i, n1, n2), work.getvector(n1, n2), t);
2746 template<
unsigned int Precision>
2748 amp::ampf<Precision>
tau,
2749 const ap::template_1d_array< amp::ampf<Precision> >&
v,
2754 ap::template_1d_array< amp::ampf<Precision> >& work)
2756 amp::ampf<Precision> t;
2761 if(
tau==0 || n1>n2 || m1>m2 )
2770 for(
i=m1;
i<=m2;
i++)
2779 for(
i=m1;
i<=m2;
i++)
2782 ap::vsub(c.getrow(
i, n1, n2),
v.getvector(1, vm), t);
2829 template<
unsigned int Precision>
2830 void rmatrixbd(ap::template_2d_array< amp::ampf<Precision> >& a,
2833 ap::template_1d_array< amp::ampf<Precision> >& tauq,
2834 ap::template_1d_array< amp::ampf<Precision> >& taup);
2835 template<
unsigned int Precision>
2836 void rmatrixbdunpackq(
const ap::template_2d_array< amp::ampf<Precision> >& qp,
2839 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
2841 ap::template_2d_array< amp::ampf<Precision> >& q);
2842 template<
unsigned int Precision>
2846 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
2847 ap::template_2d_array< amp::ampf<Precision> >& z,
2852 template<
unsigned int Precision>
2856 const ap::template_1d_array< amp::ampf<Precision> >& taup,
2858 ap::template_2d_array< amp::ampf<Precision> >& pt);
2859 template<
unsigned int Precision>
2863 const ap::template_1d_array< amp::ampf<Precision> >& taup,
2864 ap::template_2d_array< amp::ampf<Precision> >& z,
2869 template<
unsigned int Precision>
2874 ap::template_1d_array< amp::ampf<Precision> >& d,
2875 ap::template_1d_array< amp::ampf<Precision> >& e);
2876 template<
unsigned int Precision>
2877 void tobidiagonal(ap::template_2d_array< amp::ampf<Precision> >& a,
2880 ap::template_1d_array< amp::ampf<Precision> >& tauq,
2881 ap::template_1d_array< amp::ampf<Precision> >& taup);
2882 template<
unsigned int Precision>
2886 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
2888 ap::template_2d_array< amp::ampf<Precision> >& q);
2889 template<
unsigned int Precision>
2893 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
2894 ap::template_2d_array< amp::ampf<Precision> >& z,
2899 template<
unsigned int Precision>
2903 const ap::template_1d_array< amp::ampf<Precision> >& taup,
2905 ap::template_2d_array< amp::ampf<Precision> >& pt);
2906 template<
unsigned int Precision>
2910 const ap::template_1d_array< amp::ampf<Precision> >& taup,
2911 ap::template_2d_array< amp::ampf<Precision> >& z,
2916 template<
unsigned int Precision>
2921 ap::template_1d_array< amp::ampf<Precision> >& d,
2922 ap::template_1d_array< amp::ampf<Precision> >& e);
2975 template<
unsigned int Precision>
2976 void rmatrixbd(ap::template_2d_array< amp::ampf<Precision> >& a,
2979 ap::template_1d_array< amp::ampf<Precision> >& tauq,
2980 ap::template_1d_array< amp::ampf<Precision> >& taup)
2982 ap::template_1d_array< amp::ampf<Precision> > work;
2983 ap::template_1d_array< amp::ampf<Precision> > t;
2988 amp::ampf<Precision> ltau;
3005 tauq.setbounds(0, n-1);
3006 taup.setbounds(0, n-1);
3010 tauq.setbounds(0,
m-1);
3011 taup.setbounds(0,
m-1);
3019 for(
i=0;
i<=n-1;
i++)
3065 for(
i=0;
i<=
m-1;
i++)
3129 template<
unsigned int Precision>
3130 void rmatrixbdunpackq(
const ap::template_2d_array< amp::ampf<Precision> >& qp,
3133 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
3135 ap::template_2d_array< amp::ampf<Precision> >& q)
3143 if(
m==0 || n==0 || qcolumns==0 )
3151 q.setbounds(0,
m-1, 0, qcolumns-1);
3152 for(
i=0;
i<=
m-1;
i++)
3154 for(
j=0;
j<=qcolumns-1;
j++)
3203 template<
unsigned int Precision>
3207 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
3208 ap::template_2d_array< amp::ampf<Precision> >& z,
3218 ap::template_1d_array< amp::ampf<Precision> >
v;
3219 ap::template_1d_array< amp::ampf<Precision> > work;
3223 if(
m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3281 while(
i!=i2+istep );
3329 while(
i!=i2+istep );
3356 template<
unsigned int Precision>
3360 const ap::template_1d_array< amp::ampf<Precision> >& taup,
3362 ap::template_2d_array< amp::ampf<Precision> >& pt)
3370 if(
m==0 || n==0 || ptrows==0 )
3378 pt.setbounds(0, ptrows-1, 0, n-1);
3379 for(
i=0;
i<=ptrows-1;
i++)
3381 for(
j=0;
j<=n-1;
j++)
3430 template<
unsigned int Precision>
3434 const ap::template_1d_array< amp::ampf<Precision> >& taup,
3435 ap::template_2d_array< amp::ampf<Precision> >& z,
3442 ap::template_1d_array< amp::ampf<Precision> >
v;
3443 ap::template_1d_array< amp::ampf<Precision> > work;
3450 if(
m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3512 while(
i!=i2+istep );
3559 while(
i!=i2+istep );
3586 template<
unsigned int Precision>
3591 ap::template_1d_array< amp::ampf<Precision> >& d,
3592 ap::template_1d_array< amp::ampf<Precision> >& e)
3604 d.setbounds(0, n-1);
3605 e.setbounds(0, n-1);
3606 for(
i=0;
i<=n-2;
i++)
3611 d(n-1) =
b(n-1,n-1);
3615 d.setbounds(0,
m-1);
3616 e.setbounds(0,
m-1);
3617 for(
i=0;
i<=
m-2;
i++)
3622 d(
m-1) =
b(
m-1,
m-1);
3631 template<
unsigned int Precision>
3632 void tobidiagonal(ap::template_2d_array< amp::ampf<Precision> >& a,
3635 ap::template_1d_array< amp::ampf<Precision> >& tauq,
3636 ap::template_1d_array< amp::ampf<Precision> >& taup)
3638 ap::template_1d_array< amp::ampf<Precision> > work;
3639 ap::template_1d_array< amp::ampf<Precision> > t;
3643 amp::ampf<Precision> ltau;
3655 taup.setbounds(1, minmn);
3656 tauq.setbounds(1, minmn);
3762 template<
unsigned int Precision>
3766 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
3768 ap::template_2d_array< amp::ampf<Precision> >& q)
3773 ap::template_1d_array< amp::ampf<Precision> >
v;
3774 ap::template_1d_array< amp::ampf<Precision> > work;
3779 if(
m==0 || n==0 || qcolumns==0 )
3787 q.setbounds(1,
m, 1, qcolumns);
3796 for(
j=1;
j<=qcolumns;
j++)
3824 ap::vmove(
v.getvector(1, vm), qp.getcolumn(
i, ip1,
m));
3836 template<
unsigned int Precision>
3840 const ap::template_1d_array< amp::ampf<Precision> >& tauq,
3841 ap::template_2d_array< amp::ampf<Precision> >& z,
3852 ap::template_1d_array< amp::ampf<Precision> >
v;
3853 ap::template_1d_array< amp::ampf<Precision> > work;
3858 if(
m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
3917 while(
i!=i2+istep );
3955 ap::vmove(
v.getvector(1, vm), qp.getcolumn(
i, ip1,
m));
3967 while(
i!=i2+istep );
3977 template<
unsigned int Precision>
3981 const ap::template_1d_array< amp::ampf<Precision> >& taup,
3983 ap::template_2d_array< amp::ampf<Precision> >& pt)
3988 ap::template_1d_array< amp::ampf<Precision> >
v;
3989 ap::template_1d_array< amp::ampf<Precision> > work;
3994 if(
m==0 || n==0 || ptrows==0 )
4002 pt.setbounds(1, ptrows, 1, n);
4009 for(
i=1;
i<=ptrows;
i++)
4029 ap::vmove(
v.getvector(1, vm), qp.getrow(
i, ip1, n));
4051 template<
unsigned int Precision>
4055 const ap::template_1d_array< amp::ampf<Precision> >& taup,
4056 ap::template_2d_array< amp::ampf<Precision> >& z,
4064 ap::template_1d_array< amp::ampf<Precision> >
v;
4065 ap::template_1d_array< amp::ampf<Precision> > work;
4073 if(
m<=0 || n<=0 || zrows<=0 || zcolumns<=0 )
4125 ap::vmove(
v.getvector(1, vm), qp.getrow(
i, ip1, n));
4137 while(
i!=i2+istep );
4185 while(
i!=i2+istep );
4194 template<
unsigned int Precision>
4199 ap::template_1d_array< amp::ampf<Precision> >& d,
4200 ap::template_1d_array< amp::ampf<Precision> >& e)
4214 for(
i=1;
i<=n-1;
i++)
4225 for(
i=1;
i<=
m-1;
i++)
4277 template<
unsigned int Precision>
4278 void rmatrixqr(ap::template_2d_array< amp::ampf<Precision> >& a,
4281 ap::template_1d_array< amp::ampf<Precision> >&
tau);
4282 template<
unsigned int Precision>
4283 void rmatrixqrunpackq(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4286 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4288 ap::template_2d_array< amp::ampf<Precision> >& q);
4289 template<
unsigned int Precision>
4290 void rmatrixqrunpackr(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4293 ap::template_2d_array< amp::ampf<Precision> >& r);
4294 template<
unsigned int Precision>
4295 void qrdecomposition(ap::template_2d_array< amp::ampf<Precision> >& a,
4298 ap::template_1d_array< amp::ampf<Precision> >&
tau);
4299 template<
unsigned int Precision>
4300 void unpackqfromqr(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4303 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4305 ap::template_2d_array< amp::ampf<Precision> >& q);
4306 template<
unsigned int Precision>
4310 ap::template_2d_array< amp::ampf<Precision> >& q,
4311 ap::template_2d_array< amp::ampf<Precision> >& r);
4352 template<
unsigned int Precision>
4353 void rmatrixqr(ap::template_2d_array< amp::ampf<Precision> >& a,
4356 ap::template_1d_array< amp::ampf<Precision> >&
tau)
4358 ap::template_1d_array< amp::ampf<Precision> > work;
4359 ap::template_1d_array< amp::ampf<Precision> > t;
4363 amp::ampf<Precision> tmp;
4373 tau.setbounds(0, minmn-1);
4379 for(
i=0;
i<=
k-1;
i++)
4422 template<
unsigned int Precision>
4423 void rmatrixqrunpackq(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4426 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4428 ap::template_2d_array< amp::ampf<Precision> >& q)
4434 ap::template_1d_array< amp::ampf<Precision> >
v;
4435 ap::template_1d_array< amp::ampf<Precision> > work;
4439 if(
m<=0 || n<=0 || qcolumns<=0 )
4449 q.setbounds(0,
m-1, 0, qcolumns-1);
4452 for(
i=0;
i<=
m-1;
i++)
4454 for(
j=0;
j<=qcolumns-1;
j++)
4470 for(
i=
k-1;
i>=0;
i--)
4498 template<
unsigned int Precision>
4499 void rmatrixqrunpackr(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4502 ap::template_2d_array< amp::ampf<Precision> >& r)
4513 r.setbounds(0,
m-1, 0, n-1);
4514 for(
i=0;
i<=n-1;
i++)
4518 for(
i=1;
i<=
m-1;
i++)
4520 ap::vmove(r.getrow(
i, 0, n-1), r.getrow(0, 0, n-1));
4522 for(
i=0;
i<=
k-1;
i++)
4532 template<
unsigned int Precision>
4533 void qrdecomposition(ap::template_2d_array< amp::ampf<Precision> >& a,
4536 ap::template_1d_array< amp::ampf<Precision> >&
tau)
4538 ap::template_1d_array< amp::ampf<Precision> > work;
4539 ap::template_1d_array< amp::ampf<Precision> > t;
4544 amp::ampf<Precision> tmp;
4550 tau.setbounds(1, minmn);
4583 template<
unsigned int Precision>
4584 void unpackqfromqr(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4587 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4589 ap::template_2d_array< amp::ampf<Precision> >& q)
4595 ap::template_1d_array< amp::ampf<Precision> >
v;
4596 ap::template_1d_array< amp::ampf<Precision> > work;
4601 if(
m==0 || n==0 || qcolumns==0 )
4611 q.setbounds(1,
m, 1, qcolumns);
4616 for(
j=1;
j<=qcolumns;
j++)
4649 template<
unsigned int Precision>
4653 ap::template_2d_array< amp::ampf<Precision> >& q,
4654 ap::template_2d_array< amp::ampf<Precision> >& r)
4658 ap::template_1d_array< amp::ampf<Precision> >
tau;
4659 ap::template_1d_array< amp::ampf<Precision> > work;
4660 ap::template_1d_array< amp::ampf<Precision> >
v;
4670 q.setbounds(1,
m, 1,
m);
4671 r.setbounds(1,
m, 1, n);
4687 ap::vmove(r.getrow(
i, 1, n), r.getrow(1, 1, n));
4737 template<
unsigned int Precision>
4738 void rmatrixlq(ap::template_2d_array< amp::ampf<Precision> >& a,
4741 ap::template_1d_array< amp::ampf<Precision> >&
tau);
4742 template<
unsigned int Precision>
4743 void rmatrixlqunpackq(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4746 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4748 ap::template_2d_array< amp::ampf<Precision> >& q);
4749 template<
unsigned int Precision>
4750 void rmatrixlqunpackl(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4753 ap::template_2d_array< amp::ampf<Precision> >&
l);
4754 template<
unsigned int Precision>
4755 void lqdecomposition(ap::template_2d_array< amp::ampf<Precision> >& a,
4758 ap::template_1d_array< amp::ampf<Precision> >&
tau);
4759 template<
unsigned int Precision>
4760 void unpackqfromlq(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4763 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4765 ap::template_2d_array< amp::ampf<Precision> >& q);
4766 template<
unsigned int Precision>
4770 ap::template_2d_array< amp::ampf<Precision> >&
l,
4771 ap::template_2d_array< amp::ampf<Precision> >& q);
4808 template<
unsigned int Precision>
4809 void rmatrixlq(ap::template_2d_array< amp::ampf<Precision> >& a,
4812 ap::template_1d_array< amp::ampf<Precision> >&
tau)
4814 ap::template_1d_array< amp::ampf<Precision> > work;
4815 ap::template_1d_array< amp::ampf<Precision> > t;
4820 amp::ampf<Precision> tmp;
4827 tau.setbounds(0, minmn-1);
4829 for(
i=0;
i<=
k-1;
i++)
4872 template<
unsigned int Precision>
4873 void rmatrixlqunpackq(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4876 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
4878 ap::template_2d_array< amp::ampf<Precision> >& q)
4884 ap::template_1d_array< amp::ampf<Precision> >
v;
4885 ap::template_1d_array< amp::ampf<Precision> > work;
4889 if(
m<=0 || n<=0 || qrows<=0 )
4899 q.setbounds(0, qrows-1, 0, n-1);
4902 for(
i=0;
i<=qrows-1;
i++)
4904 for(
j=0;
j<=n-1;
j++)
4920 for(
i=
k-1;
i>=0;
i--)
4948 template<
unsigned int Precision>
4949 void rmatrixlqunpackl(
const ap::template_2d_array< amp::ampf<Precision> >& a,
4952 ap::template_2d_array< amp::ampf<Precision> >&
l)
4962 l.setbounds(0,
m-1, 0, n-1);
4963 for(
i=0;
i<=n-1;
i++)
4967 for(
i=1;
i<=
m-1;
i++)
4971 for(
i=0;
i<=
m-1;
i++)
4983 template<
unsigned int Precision>
4984 void lqdecomposition(ap::template_2d_array< amp::ampf<Precision> >& a,
4987 ap::template_1d_array< amp::ampf<Precision> >&
tau)
4989 ap::template_1d_array< amp::ampf<Precision> > work;
4990 ap::template_1d_array< amp::ampf<Precision> > t;
4996 amp::ampf<Precision> tmp;
5003 tau.setbounds(1, minmn);
5037 template<
unsigned int Precision>
5038 void unpackqfromlq(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5041 const ap::template_1d_array< amp::ampf<Precision> >&
tau,
5043 ap::template_2d_array< amp::ampf<Precision> >& q)
5049 ap::template_1d_array< amp::ampf<Precision> >
v;
5050 ap::template_1d_array< amp::ampf<Precision> > work;
5055 if(
m==0 || n==0 || qrows==0 )
5065 q.setbounds(1, qrows, 1, n);
5068 for(
i=1;
i<=qrows;
i++)
5103 template<
unsigned int Precision>
5107 ap::template_2d_array< amp::ampf<Precision> >&
l,
5108 ap::template_2d_array< amp::ampf<Precision> >& q)
5112 ap::template_1d_array< amp::ampf<Precision> >
tau;
5119 q.setbounds(1, n, 1, n);
5120 l.setbounds(1,
m, 1, n);
5188 template<
unsigned int Precision>
5189 amp::ampf<Precision>
vectornorm2(
const ap::template_1d_array< amp::ampf<Precision> >&
x,
5192 template<
unsigned int Precision>
5193 int vectoridxabsmax(
const ap::template_1d_array< amp::ampf<Precision> >&
x,
5196 template<
unsigned int Precision>
5197 int columnidxabsmax(
const ap::template_2d_array< amp::ampf<Precision> >&
x,
5201 template<
unsigned int Precision>
5202 int rowidxabsmax(
const ap::template_2d_array< amp::ampf<Precision> >&
x,
5206 template<
unsigned int Precision>
5207 amp::ampf<Precision>
upperhessenberg1norm(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5212 ap::template_1d_array< amp::ampf<Precision> >& work);
5213 template<
unsigned int Precision>
5214 void copymatrix(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5219 ap::template_2d_array< amp::ampf<Precision> >&
b,
5224 template<
unsigned int Precision>
5230 ap::template_1d_array< amp::ampf<Precision> >& work);
5231 template<
unsigned int Precision>
5232 void copyandtranspose(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5237 ap::template_2d_array< amp::ampf<Precision> >&
b,
5242 template<
unsigned int Precision>
5249 const ap::template_1d_array< amp::ampf<Precision> >&
x,
5252 amp::ampf<Precision>
alpha,
5253 ap::template_1d_array< amp::ampf<Precision> >&
y,
5256 amp::ampf<Precision>
beta);
5257 template<
unsigned int Precision>
5258 amp::ampf<Precision>
pythag2(amp::ampf<Precision>
x,
5259 amp::ampf<Precision>
y);
5260 template<
unsigned int Precision>
5267 const ap::template_2d_array< amp::ampf<Precision> >&
b,
5273 amp::ampf<Precision>
alpha,
5274 ap::template_2d_array< amp::ampf<Precision> >& c,
5279 amp::ampf<Precision>
beta,
5280 ap::template_1d_array< amp::ampf<Precision> >& work);
5283 template<
unsigned int Precision>
5284 amp::ampf<Precision>
vectornorm2(
const ap::template_1d_array< amp::ampf<Precision> >&
x,
5288 amp::ampf<Precision>
result;
5291 amp::ampf<Precision> absxi;
5292 amp::ampf<Precision> scl;
5293 amp::ampf<Precision> ssq;
5309 for(ix=i1; ix<=i2; ix++)
5330 template<
unsigned int Precision>
5331 int vectoridxabsmax(
const ap::template_1d_array< amp::ampf<Precision> >&
x,
5337 amp::ampf<Precision> a;
5342 for(
i=i1+1;
i<=i2;
i++)
5353 template<
unsigned int Precision>
5354 int columnidxabsmax(
const ap::template_2d_array< amp::ampf<Precision> >&
x,
5361 amp::ampf<Precision> a;
5366 for(
i=i1+1;
i<=i2;
i++)
5377 template<
unsigned int Precision>
5378 int rowidxabsmax(
const ap::template_2d_array< amp::ampf<Precision> >&
x,
5385 amp::ampf<Precision> a;
5390 for(
j=j1+1;
j<=j2;
j++)
5401 template<
unsigned int Precision>
5402 amp::ampf<Precision>
upperhessenberg1norm(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5407 ap::template_1d_array< amp::ampf<Precision> >& work)
5409 amp::ampf<Precision>
result;
5415 for(
j=j1;
j<=j2;
j++)
5419 for(
i=i1;
i<=i2;
i++)
5427 for(
j=j1;
j<=j2;
j++)
5435 template<
unsigned int Precision>
5436 void copymatrix(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5441 ap::template_2d_array< amp::ampf<Precision> >&
b,
5451 if( is1>is2 || js1>js2 )
5457 for(isrc=is1; isrc<=is2; isrc++)
5459 idst = isrc-is1+id1;
5460 ap::vmove(
b.getrow(idst, jd1, jd2), a.getrow(isrc, js1, js2));
5465 template<
unsigned int Precision>
5471 ap::template_1d_array< amp::ampf<Precision> >& work)
5480 if( i1>i2 || j1>j2 )
5485 for(
i=i1;
i<=i2-1;
i++)
5492 ap::vmove(a.getcolumn(
j, ips, i2), a.getrow(
i, jps, j2));
5498 template<
unsigned int Precision>
5499 void copyandtranspose(
const ap::template_2d_array< amp::ampf<Precision> >& a,
5504 ap::template_2d_array< amp::ampf<Precision> >&
b,
5514 if( is1>is2 || js1>js2 )
5520 for(isrc=is1; isrc<=is2; isrc++)
5522 jdst = isrc-is1+jd1;
5523 ap::vmove(
b.getcolumn(jdst, id1, id2), a.getrow(isrc, js1, js2));
5528 template<
unsigned int Precision>
5535 const ap::template_1d_array< amp::ampf<Precision> >&
x,
5538 amp::ampf<Precision>
alpha,
5539 ap::template_1d_array< amp::ampf<Precision> >&
y,
5542 amp::ampf<Precision>
beta)
5545 amp::ampf<Precision>
v;
5554 if( i1>i2 || j1>j2 )
5566 for(
i=iy1;
i<=iy2;
i++)
5579 for(
i=i1;
i<=i2;
i++)
5591 if( i1>i2 || j1>j2 )
5603 for(
i=iy1;
i<=iy2;
i++)
5616 for(
i=i1;
i<=i2;
i++)
5619 ap::vadd(
y.getvector(iy1, iy2), a.getrow(
i, j1, j2),
v);
5625 template<
unsigned int Precision>
5626 amp::ampf<Precision>
pythag2(amp::ampf<Precision>
x,
5627 amp::ampf<Precision>
y)
5629 amp::ampf<Precision>
result;
5630 amp::ampf<Precision>
w;
5631 amp::ampf<Precision> xabs;
5632 amp::ampf<Precision> yabs;
5633 amp::ampf<Precision> z;
5652 template<
unsigned int Precision>
5659 const ap::template_2d_array< amp::ampf<Precision> >&
b,
5665 amp::ampf<Precision>
alpha,
5666 ap::template_2d_array< amp::ampf<Precision> >& c,
5671 amp::ampf<Precision>
beta,
5672 ap::template_1d_array< amp::ampf<Precision> >& work)
5685 amp::ampf<Precision>
v;
5713 if( arows<=0 || acols<=0 || brows<=0 || bcols<=0 )
5734 for(
i=ci1;
i<=ci2;
i++)
5736 for(
j=cj1;
j<=cj2;
j++)
5744 for(
i=ci1;
i<=ci2;
i++)
5753 if( !transa && !transb )
5755 for(
l=ai1;
l<=ai2;
l++)
5757 for(r=bi1; r<=bi2; r++)
5761 ap::vadd(c.getrow(
k, cj1, cj2),
b.getrow(r, bj1, bj2),
v);
5770 if( !transa && transb )
5772 if( arows*acols<brows*bcols )
5774 for(r=bi1; r<=bi2; r++)
5776 for(
l=ai1;
l<=ai2;
l++)
5779 c(ci1+
l-ai1,cj1+r-bi1) = c(ci1+
l-ai1,cj1+r-bi1)+
alpha*
v;
5786 for(
l=ai1;
l<=ai2;
l++)
5788 for(r=bi1; r<=bi2; r++)
5791 c(ci1+
l-ai1,cj1+r-bi1) = c(ci1+
l-ai1,cj1+r-bi1)+
alpha*
v;
5801 if( transa && !transb )
5803 for(
l=aj1;
l<=aj2;
l++)
5805 for(r=bi1; r<=bi2; r++)
5809 ap::vadd(c.getrow(
k, cj1, cj2),
b.getrow(r, bj1, bj2),
v);
5818 if( transa && transb )
5820 if( arows*acols<brows*bcols )
5822 for(r=bi1; r<=bi2; r++)
5824 for(
i=1;
i<=crows;
i++)
5826 work(
i) = amp::ampf<Precision>(
"0.0");
5828 for(
l=ai1;
l<=ai2;
l++)
5840 for(
l=aj1;
l<=aj2;
l++)
5844 for(r=bi1; r<=bi2; r++)
5847 c(ci1+
l-aj1,cj1+r-bi1) = c(ci1+
l-aj1,cj1+r-bi1)+
alpha*
v;
5898 template<
unsigned int Precision>
5904 const ap::template_1d_array< amp::ampf<Precision> >& c,
5905 const ap::template_1d_array< amp::ampf<Precision> >&
s,
5906 ap::template_2d_array< amp::ampf<Precision> >& a,
5907 ap::template_1d_array< amp::ampf<Precision> >& work);
5908 template<
unsigned int Precision>
5914 const ap::template_1d_array< amp::ampf<Precision> >& c,
5915 const ap::template_1d_array< amp::ampf<Precision> >&
s,
5916 ap::template_2d_array< amp::ampf<Precision> >& a,
5917 ap::template_1d_array< amp::ampf<Precision> >& work);
5918 template<
unsigned int Precision>
5920 amp::ampf<Precision>
g,
5921 amp::ampf<Precision>& cs,
5922 amp::ampf<Precision>& sn,
5923 amp::ampf<Precision>& r);
5951 template<
unsigned int Precision>
5957 const ap::template_1d_array< amp::ampf<Precision> >& c,
5958 const ap::template_1d_array< amp::ampf<Precision> >&
s,
5959 ap::template_2d_array< amp::ampf<Precision> >& a,
5960 ap::template_1d_array< amp::ampf<Precision> >& work)
5964 amp::ampf<Precision> ctemp;
5965 amp::ampf<Precision> stemp;
5966 amp::ampf<Precision> temp;
5969 if( m1>m2 || n1>n2 )
5985 for(
j=m1;
j<=m2-1;
j++)
5989 if( ctemp!=1 || stemp!=0 )
5995 ap::vadd(a.getrow(
j, n1, n2), a.getrow(jp1, n1, n2), stemp);
6006 for(
j=m1;
j<=m2-1;
j++)
6010 if( ctemp!=1 || stemp!=0 )
6013 a(
j+1,n1) = ctemp*temp-stemp*a(
j,n1);
6014 a(
j,n1) = stemp*temp+ctemp*a(
j,n1);
6027 for(
j=m2-1;
j>=m1;
j--)
6031 if( ctemp!=1 || stemp!=0 )
6037 ap::vadd(a.getrow(
j, n1, n2), a.getrow(jp1, n1, n2), stemp);
6048 for(
j=m2-1;
j>=m1;
j--)
6052 if( ctemp!=1 || stemp!=0 )
6055 a(
j+1,n1) = ctemp*temp-stemp*a(
j,n1);
6056 a(
j,n1) = stemp*temp+ctemp*a(
j,n1);
6089 template<
unsigned int Precision>
6095 const ap::template_1d_array< amp::ampf<Precision> >& c,
6096 const ap::template_1d_array< amp::ampf<Precision> >&
s,
6097 ap::template_2d_array< amp::ampf<Precision> >& a,
6098 ap::template_1d_array< amp::ampf<Precision> >& work)
6102 amp::ampf<Precision> ctemp;
6103 amp::ampf<Precision> stemp;
6104 amp::ampf<Precision> temp;
6119 for(
j=n1;
j<=n2-1;
j++)
6123 if( ctemp!=1 || stemp!=0 )
6128 ap::vmul(a.getcolumn(
j, m1, m2), ctemp);
6129 ap::vadd(a.getcolumn(
j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
6140 for(
j=n1;
j<=n2-1;
j++)
6144 if( ctemp!=1 || stemp!=0 )
6147 a(m1,
j+1) = ctemp*temp-stemp*a(m1,
j);
6148 a(m1,
j) = stemp*temp+ctemp*a(m1,
j);
6161 for(
j=n2-1;
j>=n1;
j--)
6165 if( ctemp!=1 || stemp!=0 )
6170 ap::vmul(a.getcolumn(
j, m1, m2), ctemp);
6171 ap::vadd(a.getcolumn(
j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
6182 for(
j=n2-1;
j>=n1;
j--)
6186 if( ctemp!=1 || stemp!=0 )
6189 a(m1,
j+1) = ctemp*temp-stemp*a(m1,
j);
6190 a(m1,
j) = stemp*temp+ctemp*a(m1,
j);
6206 template<
unsigned int Precision>
6208 amp::ampf<Precision>
g,
6209 amp::ampf<Precision>& cs,
6210 amp::ampf<Precision>& sn,
6211 amp::ampf<Precision>& r)
6213 amp::ampf<Precision> f1;
6214 amp::ampf<Precision> g1;
6291 template<
unsigned int Precision>
6292 bool rmatrixbdsvd(ap::template_1d_array< amp::ampf<Precision> >& d,
6293 ap::template_1d_array< amp::ampf<Precision> > e,
6296 bool isfractionalaccuracyrequired,
6297 ap::template_2d_array< amp::ampf<Precision> >& u,
6299 ap::template_2d_array< amp::ampf<Precision> >& c,
6301 ap::template_2d_array< amp::ampf<Precision> >& vt,
6303 template<
unsigned int Precision>
6305 ap::template_1d_array< amp::ampf<Precision> > e,
6308 bool isfractionalaccuracyrequired,
6309 ap::template_2d_array< amp::ampf<Precision> >& u,
6311 ap::template_2d_array< amp::ampf<Precision> >& c,
6313 ap::template_2d_array< amp::ampf<Precision> >& vt,
6315 template<
unsigned int Precision>
6317 ap::template_1d_array< amp::ampf<Precision> > e,
6320 bool isfractionalaccuracyrequired,
6321 ap::template_2d_array< amp::ampf<Precision> >& u,
6324 ap::template_2d_array< amp::ampf<Precision> >& c,
6327 ap::template_2d_array< amp::ampf<Precision> >& vt,
6330 template<
unsigned int Precision>
6331 amp::ampf<Precision>
extsignbdsqr(amp::ampf<Precision> a,
6332 amp::ampf<Precision>
b);
6333 template<
unsigned int Precision>
6334 void svd2x2(amp::ampf<Precision>
f,
6335 amp::ampf<Precision>
g,
6336 amp::ampf<Precision>
h,
6337 amp::ampf<Precision>& ssmin,
6338 amp::ampf<Precision>& ssmax);
6339 template<
unsigned int Precision>
6340 void svdv2x2(amp::ampf<Precision>
f,
6341 amp::ampf<Precision>
g,
6342 amp::ampf<Precision>
h,
6343 amp::ampf<Precision>& ssmin,
6344 amp::ampf<Precision>& ssmax,
6345 amp::ampf<Precision>& snr,
6346 amp::ampf<Precision>& csr,
6347 amp::ampf<Precision>& snl,
6348 amp::ampf<Precision>& csl);
6429 template<
unsigned int Precision>
6430 bool rmatrixbdsvd(ap::template_1d_array< amp::ampf<Precision> >& d,
6431 ap::template_1d_array< amp::ampf<Precision> > e,
6434 bool isfractionalaccuracyrequired,
6435 ap::template_2d_array< amp::ampf<Precision> >& u,
6437 ap::template_2d_array< amp::ampf<Precision> >& c,
6439 ap::template_2d_array< amp::ampf<Precision> >& vt,
6443 ap::template_1d_array< amp::ampf<Precision> > d1;
6444 ap::template_1d_array< amp::ampf<Precision> > e1;
6454 result =
bidiagonalsvddecompositioninternal<Precision>(d1, e1, n, isupper, isfractionalaccuracyrequired, u, 0, nru, c, 0, ncc, vt, 0, ncvt);
6472 template<
unsigned int Precision>
6474 ap::template_1d_array< amp::ampf<Precision> > e,
6477 bool isfractionalaccuracyrequired,
6478 ap::template_2d_array< amp::ampf<Precision> >& u,
6480 ap::template_2d_array< amp::ampf<Precision> >& c,
6482 ap::template_2d_array< amp::ampf<Precision> >& vt,
6488 result =
bidiagonalsvddecompositioninternal<Precision>(d, e, n, isupper, isfractionalaccuracyrequired, u, 1, nru, c, 1, ncc, vt, 1, ncvt);
6496 template<
unsigned int Precision>
6498 ap::template_1d_array< amp::ampf<Precision> > e,
6501 bool isfractionalaccuracyrequired,
6502 ap::template_2d_array< amp::ampf<Precision> >& u,
6505 ap::template_2d_array< amp::ampf<Precision> >& c,
6508 ap::template_2d_array< amp::ampf<Precision> >& vt,
6524 amp::ampf<Precision> abse;
6525 amp::ampf<Precision> abss;
6526 amp::ampf<Precision> cosl;
6527 amp::ampf<Precision> cosr;
6528 amp::ampf<Precision> cs;
6529 amp::ampf<Precision> eps;
6530 amp::ampf<Precision>
f;
6531 amp::ampf<Precision>
g;
6532 amp::ampf<Precision>
h;
6533 amp::ampf<Precision>
mu;
6534 amp::ampf<Precision> oldcs;
6535 amp::ampf<Precision> oldsn;
6536 amp::ampf<Precision> r;
6537 amp::ampf<Precision> shift;
6538 amp::ampf<Precision> sigmn;
6539 amp::ampf<Precision> sigmx;
6540 amp::ampf<Precision> sinl;
6541 amp::ampf<Precision> sinr;
6542 amp::ampf<Precision> sll;
6543 amp::ampf<Precision> smax;
6544 amp::ampf<Precision> smin;
6545 amp::ampf<Precision> sminl;
6546 amp::ampf<Precision> sminlo;
6547 amp::ampf<Precision> sminoa;
6548 amp::ampf<Precision> sn;
6549 amp::ampf<Precision> thresh;
6550 amp::ampf<Precision> tol;
6551 amp::ampf<Precision> tolmul;
6552 amp::ampf<Precision> unfl;
6553 ap::template_1d_array< amp::ampf<Precision> > work0;
6554 ap::template_1d_array< amp::ampf<Precision> > work1;
6555 ap::template_1d_array< amp::ampf<Precision> > work2;
6556 ap::template_1d_array< amp::ampf<Precision> > work3;
6558 bool matrixsplitflag;
6560 ap::template_1d_array< amp::ampf<Precision> > utemp;
6561 ap::template_1d_array< amp::ampf<Precision> > vttemp;
6562 ap::template_1d_array< amp::ampf<Precision> > ctemp;
6563 ap::template_1d_array< amp::ampf<Precision> > etemp;
6565 amp::ampf<Precision> tmp;
6586 ap::vmul(vt.getrow(vstart, vstart, vstart+ncvt-1), -1);
6612 for(
i=1;
i<=n-1;
i++)
6617 for(
i=1;
i<=n-1;
i++)
6636 for(
i=1;
i<=n-1;
i++)
6666 if( !isfractionalaccuracyrequired )
6679 for(
i=1;
i<=n-1;
i++)
6759 matrixsplitflag =
false;
6760 for(lll=1; lll<=
m-1; lll++)
6765 if( tol<0 && abss<=thresh )
6771 matrixsplitflag =
true;
6777 if( !matrixsplitflag )
6820 mm1 =
m-1+(vstart-1);
6823 ap::vmul(vt.getrow(mm0, vstart, vend), cosr);
6824 ap::vsub(vt.getrow(mm0, vstart, vend), vt.getrow(mm1, vstart, vend), sinr);
6833 ap::vmul(u.getcolumn(mm0, ustart, uend), cosl);
6834 ap::vsub(u.getcolumn(mm0, ustart, uend), u.getcolumn(mm1, ustart, uend), sinl);
6843 ap::vmul(c.getrow(mm0, cstart, cend), cosl);
6844 ap::vsub(c.getrow(mm0, cstart, cend), c.getrow(mm1, cstart, cend), sinl);
6869 if( ll!=oldll ||
m!=oldm || bchangedir )
6914 for(lll=ll; lll<=
m-1; lll++)
6954 for(lll=
m-1; lll>=ll; lll--)
7035 for(
i=ll;
i<=
m-1;
i++)
7046 work2(
i-ll+1) = oldcs;
7047 work3(
i-ll+1) = oldsn;
7086 for(
i=
m;
i>=ll+1;
i--)
7097 work2(
i-ll) = oldcs;
7098 work3(
i-ll) = -oldsn;
7144 for(
i=ll;
i<=
m-1;
i++)
7151 f = cosr*d(
i)+sinr*e(
i);
7152 e(
i) = cosr*e(
i)-sinr*d(
i);
7154 d(
i+1) = cosr*d(
i+1);
7157 f = cosl*e(
i)+sinl*d(
i+1);
7158 d(
i+1) = cosl*d(
i+1)-sinl*e(
i);
7162 e(
i+1) = cosl*e(
i+1);
7164 work0(
i-ll+1) = cosr;
7165 work1(
i-ll+1) = sinr;
7166 work2(
i-ll+1) = cosl;
7167 work3(
i-ll+1) = sinl;
7204 for(
i=
m;
i>=ll+1;
i--)
7211 f = cosr*d(
i)+sinr*e(
i-1);
7212 e(
i-1) = cosr*e(
i-1)-sinr*d(
i);
7214 d(
i-1) = cosr*d(
i-1);
7217 f = cosl*e(
i-1)+sinl*d(
i-1);
7218 d(
i-1) = cosl*d(
i-1)-sinl*e(
i-1);
7222 e(
i-2) = cosl*e(
i-2);
7225 work1(
i-ll) = -sinr;
7227 work3(
i-ll) = -sinl;
7277 ap::vmul(vt.getrow(
i+vstart-1, vstart, vend), -1);
7286 for(
i=1;
i<=n-1;
i++)
7294 for(
j=2;
j<=n+1-
i;
j++)
7314 ap::vmove(vt.getrow(isub+vstart-1, vstart, vend), vt.getrow(
j+vstart-1, vstart, vend));
7321 ap::vmove(u.getcolumn(isub+ustart-1, ustart, uend), u.getcolumn(
j+ustart-1, ustart, uend));
7328 ap::vmove(c.getrow(isub+cstart-1, cstart, cend), c.getrow(
j+cstart-1, cstart, cend));
7337 template<
unsigned int Precision>
7338 amp::ampf<Precision>
extsignbdsqr(amp::ampf<Precision> a,
7339 amp::ampf<Precision>
b)
7341 amp::ampf<Precision>
result;
7356 template<
unsigned int Precision>
7357 void svd2x2(amp::ampf<Precision>
f,
7358 amp::ampf<Precision>
g,
7359 amp::ampf<Precision>
h,
7360 amp::ampf<Precision>& ssmin,
7361 amp::ampf<Precision>& ssmax)
7363 amp::ampf<Precision> aas;
7364 amp::ampf<Precision> at;
7365 amp::ampf<Precision> au;
7366 amp::ampf<Precision> c;
7367 amp::ampf<Precision>
fa;
7368 amp::ampf<Precision> fhmn;
7369 amp::ampf<Precision> fhmx;
7370 amp::ampf<Precision> ga;
7371 amp::ampf<Precision> ha;
7396 at = (fhmx-fhmn)/fhmx;
7413 ssmin = fhmn*fhmx/ga;
7419 at = (fhmx-fhmn)/fhmx;
7422 ssmin = ssmin+ssmin;
7430 template<
unsigned int Precision>
7431 void svdv2x2(amp::ampf<Precision>
f,
7432 amp::ampf<Precision>
g,
7433 amp::ampf<Precision>
h,
7434 amp::ampf<Precision>& ssmin,
7435 amp::ampf<Precision>& ssmax,
7436 amp::ampf<Precision>& snr,
7437 amp::ampf<Precision>& csr,
7438 amp::ampf<Precision>& snl,
7439 amp::ampf<Precision>& csl)
7444 amp::ampf<Precision> a;
7445 amp::ampf<Precision> clt;
7446 amp::ampf<Precision> crt;
7447 amp::ampf<Precision> d;
7448 amp::ampf<Precision>
fa;
7449 amp::ampf<Precision> ft;
7450 amp::ampf<Precision> ga;
7451 amp::ampf<Precision> gt;
7452 amp::ampf<Precision> ha;
7453 amp::ampf<Precision> ht;
7454 amp::ampf<Precision>
l;
7455 amp::ampf<Precision>
m;
7456 amp::ampf<Precision> mm;
7457 amp::ampf<Precision> r;
7458 amp::ampf<Precision>
s;
7459 amp::ampf<Precision> slt;
7460 amp::ampf<Precision> srt;
7461 amp::ampf<Precision> t;
7462 amp::ampf<Precision> temp;
7463 amp::ampf<Precision> tsign;
7464 amp::ampf<Precision> tt;
7465 amp::ampf<Precision>
v;
7568 a = amp::ampf<Precision>(
"0.5")*(
s+r);
7588 t = (
m/(
s+t)+
m/(r+
l))*(1+a);
7593 clt = (crt+srt*
m)/a;
7687 template<
unsigned int Precision>
7688 bool rmatrixsvd(ap::template_2d_array< amp::ampf<Precision> > a,
7693 int additionalmemory,
7694 ap::template_1d_array< amp::ampf<Precision> >&
w,
7695 ap::template_2d_array< amp::ampf<Precision> >& u,
7696 ap::template_2d_array< amp::ampf<Precision> >& vt);
7697 template<
unsigned int Precision>
7703 int additionalmemory,
7704 ap::template_1d_array< amp::ampf<Precision> >&
w,
7705 ap::template_2d_array< amp::ampf<Precision> >& u,
7706 ap::template_2d_array< amp::ampf<Precision> >& vt);
7759 template<
unsigned int Precision>
7760 bool rmatrixsvd(ap::template_2d_array< amp::ampf<Precision> > a,
7765 int additionalmemory,
7766 ap::template_1d_array< amp::ampf<Precision> >&
w,
7767 ap::template_2d_array< amp::ampf<Precision> >& u,
7768 ap::template_2d_array< amp::ampf<Precision> >& vt)
7771 ap::template_1d_array< amp::ampf<Precision> > tauq;
7772 ap::template_1d_array< amp::ampf<Precision> > taup;
7773 ap::template_1d_array< amp::ampf<Precision> >
tau;
7774 ap::template_1d_array< amp::ampf<Precision> > e;
7775 ap::template_1d_array< amp::ampf<Precision> > work;
7776 ap::template_2d_array< amp::ampf<Precision> > t2;
7786 amp::ampf<Precision> sm;
7802 w.setbounds(1, minmn);
7809 u.setbounds(0, nru-1, 0, ncu-1);
7815 u.setbounds(0, nru-1, 0, ncu-1);
7823 vt.setbounds(0, nrvt-1, 0, ncvt-1);
7829 vt.setbounds(0, nrvt-1, 0, ncvt-1);
7836 if(
m>amp::ampf<Precision>(
"1.6")*n )
7845 for(
i=0;
i<=n-1;
i++)
7847 for(
j=0;
j<=
i-1;
j++)
7855 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u, 0, a, 0, vt, ncvt);
7866 for(
i=0;
i<=n-1;
i++)
7868 for(
j=0;
j<=
i-1;
j++)
7876 if( additionalmemory<1 )
7883 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u,
m, a, 0, vt, ncvt);
7895 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e, n, isupper,
false, u, 0, t2, n, vt, ncvt);
7896 blas::matrixmatrixmultiply<Precision>(a, 0,
m-1, 0, n-1,
false, t2, 0, n-1, 0, n-1,
true, amp::ampf<Precision>(
"1.0"), u, 0,
m-1, 0, n-1, amp::ampf<Precision>(
"0.0"), work);
7906 if( n>amp::ampf<Precision>(
"1.6")*
m )
7915 for(
i=0;
i<=
m-1;
i++)
7917 for(
j=
i+1;
j<=
m-1;
j++)
7927 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, 0);
7939 for(
i=0;
i<=
m-1;
i++)
7941 for(
j=
i+1;
j<=
m-1;
j++)
7951 if( additionalmemory<1 )
7958 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, n);
7967 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, t2,
m);
7969 blas::matrixmatrixmultiply<Precision>(t2, 0,
m-1, 0,
m-1,
false, a, 0,
m-1, 0, n-1,
false, amp::ampf<Precision>(
"1.0"), vt, 0,
m-1, 0, n-1, amp::ampf<Precision>(
"0.0"), work);
7988 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, a, 0, u, nru, vt, ncvt);
8000 if( additionalmemory<2 || uneeded==0 )
8006 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, u, nru, a, 0, vt, ncvt);
8016 result =
bdsvd::rmatrixbdsvd<Precision>(
w, e, minmn, isupper,
false, u, 0, t2,
m, vt, ncvt);
8027 template<
unsigned int Precision>
8033 int additionalmemory,
8034 ap::template_1d_array< amp::ampf<Precision> >&
w,
8035 ap::template_2d_array< amp::ampf<Precision> >& u,
8036 ap::template_2d_array< amp::ampf<Precision> >& vt)
8039 ap::template_1d_array< amp::ampf<Precision> > tauq;
8040 ap::template_1d_array< amp::ampf<Precision> > taup;
8041 ap::template_1d_array< amp::ampf<Precision> >
tau;
8042 ap::template_1d_array< amp::ampf<Precision> > e;
8043 ap::template_1d_array< amp::ampf<Precision> > work;
8044 ap::template_2d_array< amp::ampf<Precision> > t2;
8054 amp::ampf<Precision> sm;
8070 w.setbounds(1, minmn);
8077 u.setbounds(1, nru, 1, ncu);
8083 u.setbounds(1, nru, 1, ncu);
8091 vt.setbounds(1, nrvt, 1, ncvt);
8097 vt.setbounds(1, nrvt, 1, ncvt);
8104 if(
m>amp::ampf<Precision>(
"1.6")*n )
8115 for(
j=1;
j<=
i-1;
j++)
8123 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u, 0, a, 0, vt, ncvt);
8136 for(
j=1;
j<=
i-1;
j++)
8144 if( additionalmemory<1 )
8151 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u,
m, a, 0, vt, ncvt);
8163 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, n, isupper,
false, u, 0, t2, n, vt, ncvt);
8164 blas::matrixmatrixmultiply<Precision>(a, 1,
m, 1, n,
false, t2, 1, n, 1, n,
true, amp::ampf<Precision>(
"1.0"), u, 1,
m, 1, n, amp::ampf<Precision>(
"0.0"), work);
8174 if( n>amp::ampf<Precision>(
"1.6")*
m )
8183 for(
i=1;
i<=
m-1;
i++)
8185 for(
j=
i+1;
j<=
m;
j++)
8195 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, 0);
8207 for(
i=1;
i<=
m-1;
i++)
8209 for(
j=
i+1;
j<=
m;
j++)
8219 if( additionalmemory<1 )
8226 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, vt, n);
8235 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e,
m, isupper,
false, a, 0, u, nru, t2,
m);
8237 blas::matrixmatrixmultiply<Precision>(t2, 1,
m, 1,
m,
false, a, 1,
m, 1, n,
false, amp::ampf<Precision>(
"1.0"), vt, 1,
m, 1, n, amp::ampf<Precision>(
"0.0"), work);
8256 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, a, 0, u, nru, vt, ncvt);
8268 if( additionalmemory<2 || uneeded==0 )
8274 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, u, nru, a, 0, vt, ncvt);
8284 result =
bdsvd::bidiagonalsvddecomposition<Precision>(
w, e, minmn, isupper,
false, u, 0, t2,
m, vt, ncvt);
#define __AMP_BINARY_OPF(type)
#define __AMP_BINARY_OPI(type)
void tau(int **points, int sizePoints, int k)
ampf & operator/=(const T &v)
ampf & operator-=(const T &v)
static const ampf getAlgoPascalEpsilon()
static const ampf getMaxNumber()
void InitializeAsSLong(signed long v)
void InitializeAsULong(unsigned long v)
static const ampf getAlgoPascalMaxNumber()
static const ampf getUlp256()
bool isNegativeNumber() const
void InitializeAsDouble(long double v)
static const ampf getAlgoPascalMinNumber()
static const ampf getMinNumber()
static const ampf getUlp()
ampf & operator+=(const T &v)
static const ampf getAlgoPascalMaxNumber()
static const ampf getRandom()
bool isPositiveNumber() const
void InitializeAsString(const char *s)
static const ampf getUlp512()
mpfr_srcptr getReadPtr() const
static const ampf getUlp512()
static const ampf getRandom()
static const ampf getAlgoPascalEpsilon()
ampf & operator=(long double v)
static const ampf getUlpOf(const ampf &x)
static const ampf getUlp()
std::string toHex() const
ampf & operator*=(const T &v)
static const ampf getMinNumber()
bool isFiniteNumber() const
ampf(const ampf< Precision2 > &r)
std::string toDec() const
ampf(const std::string &s)
static const ampf getMaxNumber()
static const ampf getUlp256()
static const ampf getAlgoPascalMinNumber()
campf & operator=(long double v)
campf(const ampf< Precision > &_x)
campf(const ampf< Precision > &_x, const ampf< Precision > &_y)
campf(const campf< Prec2 > &z)
void initialize(int Precision)
mpfr_reference(const mpfr_reference &r)
mpfr_reference & operator=(const mpfr_reference &r)
mpfr_srcptr getReadPtr() const
static gmp_randstate_t * getRandState()
static gmp_randstate_t * getRandState()
static void deleteMpfr(mpfr_record *ref)
static mpfr_record * newMpfr(unsigned int Precision)
static void deleteMpfr(mpfr_record *ref)
static mpfr_record_ptr & getList(unsigned int Precision)
static mpfr_record * newMpfr(unsigned int Precision)
static void make_assertion(bool bClause)
complex & operator+=(const double &v)
complex & operator+=(const complex &z)
complex(const double &_x)
complex & operator-=(const complex &z)
complex & operator=(const double &v)
complex & operator-=(const double &v)
complex & operator*=(const complex &z)
complex & operator*=(const double &v)
complex & operator/=(const double &v)
complex & operator/=(const complex &z)
complex(const double &_x, const double &_y)
complex(const complex &z)
const T * GetData() const
const_raw_vector(const T *Data, int Length, int Step)
raw_vector(T *Data, int Length, int Step)
raw_vector< T > getvector(int iStart, int iEnd)
int gethighbound(int iBoundNum=0) const
bool wrongIdx(int i) const
const_raw_vector< T > getvector(int iStart, int iEnd) const
const T & operator()(int i) const
int getlowbound(int iBoundNum=0) const
void setcontent(int iLow, int iHigh, const T *pContent)
template_1d_array(const template_1d_array &rhs)
void setbounds(int iLow, int iHigh)
const template_1d_array & operator=(const template_1d_array &rhs)
const T * getcontent() const
void setcontent(int iLow1, int iHigh1, int iLow2, int iHigh2, const T *pContent)
template_2d_array(const template_2d_array &rhs)
const T & operator()(int i1, int i2) const
void setbounds(int iLow1, int iHigh1, int iLow2, int iHigh2)
bool wrongRow(int i) const
bool wrongColumn(int j) const
const T * getcontent() const
int getlowbound(int iBoundNum) const
const_raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd) const
const template_2d_array & operator=(const template_2d_array &rhs)
raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd)
raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd)
int gethighbound(int iBoundNum) const
const_raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd) const
T & operator()(int i1, int i2)
static BOOLEAN fa(leftv res, leftv args)
const CanonicalForm int s
const CanonicalForm int const CFList const Variable & y
const Variable & v
< [in] a sqrfree bivariate poly
void WerrorS(const char *s)
static matrix mu(matrix A, const ring R)
campf< Precision > & operator-=(campf< Precision > &lhs, const campf< Precision > &rhs)
campf< Precision > & operator*=(campf< Precision > &lhs, const campf< Precision > &rhs)
const ampf< Precision > abs(const ampf< Precision > &x)
const ampf< Precision > cos(const ampf< Precision > &x)
const ampf< Precision > halfpi()
campf< Precision > & operator/=(campf< Precision > &lhs, const campf< Precision > &rhs)
const signed long floor(const ampf< Precision > &x)
const bool operator>(const ampf< Precision > &op1, const ampf< Precision > &op2)
const ampf< Precision > operator-(const ampf< Precision > &op1)
void vAdd(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
const signed long ceil(const ampf< Precision > &x)
const ampf< Precision > abscomplex(const campf< Precision > &z)
const ampf< Precision > operator/(const ampf< Precision > &op1, const ampf< Precision > &op2)
const ampf< Precision > cosh(const ampf< Precision > &x)
const ampf< Precision > ldexp2(const ampf< Precision > &x, mp_exp_t exponent)
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
const ampf< Precision > exp(const ampf< Precision > &x)
const ampf< Precision > sqrt(const ampf< Precision > &x)
const campf< Precision > conj(const campf< Precision > &z)
const ampf< Precision > twopi()
const bool operator>=(const ampf< Precision > &op1, const ampf< Precision > &op2)
void vMoveNeg(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
const ampf< Precision > atan2(const ampf< Precision > &y, const ampf< Precision > &x)
const campf< Precision > csqr(const campf< Precision > &z)
const ampf< Precision > log2(const ampf< Precision > &x)
const bool operator<=(const ampf< Precision > &op1, const ampf< Precision > &op2)
const ampf< Precision > operator+(const ampf< Precision > &op1)
const bool operator!=(const ampf< Precision > &op1, const ampf< Precision > &op2)
void vMove(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
const ampf< Precision > log10(const ampf< Precision > &x)
const ampf< Precision > asin(const ampf< Precision > &x)
const bool operator<(const ampf< Precision > &op1, const ampf< Precision > &op2)
const ampf< Precision > pow(const ampf< Precision > &x, const ampf< Precision > &y)
campf< Precision > & operator+=(campf< Precision > &lhs, const campf< Precision > &rhs)
const ampf< Precision > tanh(const ampf< Precision > &x)
const ampf< Precision > frexp2(const ampf< Precision > &x, mp_exp_t *exponent)
void vSub(ap::raw_vector< ampf< Precision > > vDst, ap::const_raw_vector< ampf< Precision > > vSrc)
const ampf< Precision > atan(const ampf< Precision > &x)
void vMul(ap::raw_vector< ampf< Precision > > vDst, T2 alpha)
const int sign(const ampf< Precision > &x)
const ampf< Precision > sin(const ampf< Precision > &x)
mpfr_record * mpfr_record_ptr
const ampf< Precision > sqr(const ampf< Precision > &x)
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
const ampf< Precision > sinh(const ampf< Precision > &x)
const ampf< Precision > operator*(const ampf< Precision > &op1, const ampf< Precision > &op2)
const signed long trunc(const ampf< Precision > &x)
const ampf< Precision > frac(const ampf< Precision > &x)
const bool operator==(const ampf< Precision > &op1, const ampf< Precision > &op2)
const ampf< Precision > acos(const ampf< Precision > &x)
const ampf< Precision > tan(const ampf< Precision > &x)
const signed long round(const ampf< Precision > &x)
const ampf< Precision > log(const ampf< Precision > &x)
int maxint(int m1, int m2)
template_2d_array< double > real_2d_array
template_1d_array< bool > boolean_1d_array
template_2d_array< bool > boolean_2d_array
const complex operator*(const complex &lhs, const complex &rhs)
const complex csqr(const complex &z)
T vdotproduct(const_raw_vector< T > v1, const_raw_vector< T > v2)
const double machineepsilon
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
template_2d_array< int > integer_2d_array
const double minrealnumber
template_2d_array< complex > complex_2d_array
void vmoveneg(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
int randominteger(int maxv)
double maxreal(double m1, double m2)
const double abscomplex(const complex &z)
template_1d_array< complex > complex_1d_array
const bool operator!=(const complex &lhs, const complex &rhs)
double minreal(double m1, double m2)
int minint(int m1, int m2)
void vmul(raw_vector< T > vdst, T2 alpha)
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
const complex conj(const complex &z)
const complex operator+(const complex &lhs)
const complex operator/(const complex &lhs, const complex &rhs)
template_1d_array< int > integer_1d_array
const double maxrealnumber
const bool operator==(const complex &lhs, const complex &rhs)
const complex operator-(const complex &lhs)
template_1d_array< double > real_1d_array
bool bidiagonalsvddecomposition(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int ncvt)
amp::ampf< Precision > extsignbdsqr(amp::ampf< Precision > a, amp::ampf< Precision > b)
bool bidiagonalsvddecompositioninternal(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int ustart, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int cstart, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int vstart, int ncvt)
bool rmatrixbdsvd(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int ncvt)
void svdv2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax, amp::ampf< Precision > &snr, amp::ampf< Precision > &csr, amp::ampf< Precision > &snl, amp::ampf< Precision > &csl)
void svd2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax)
void unpackptfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
void multiplybyqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
void rmatrixbdunpackpt(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, int ptrows, ap::template_2d_array< amp::ampf< Precision > > &pt)
void unpackqfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
void rmatrixbdunpackq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
void rmatrixbdmultiplybyq(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
void unpackdiagonalsfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
void rmatrixbdunpackdiagonals(const ap::template_2d_array< amp::ampf< Precision > > &b, int m, int n, bool &isupper, ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > &e)
void tobidiagonal(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
void rmatrixbd(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tauq, ap::template_1d_array< amp::ampf< Precision > > &taup)
void rmatrixbdmultiplybyp(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
void multiplybypfrombidiagonal(const ap::template_2d_array< amp::ampf< Precision > > &qp, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &taup, ap::template_2d_array< amp::ampf< Precision > > &z, int zrows, int zcolumns, bool fromtheright, bool dotranspose)
int columnidxabsmax(const ap::template_2d_array< amp::ampf< Precision > > &x, int i1, int i2, int j)
void matrixmatrixmultiply(const ap::template_2d_array< amp::ampf< Precision > > &a, int ai1, int ai2, int aj1, int aj2, bool transa, const ap::template_2d_array< amp::ampf< Precision > > &b, int bi1, int bi2, int bj1, int bj2, bool transb, amp::ampf< Precision > alpha, ap::template_2d_array< amp::ampf< Precision > > &c, int ci1, int ci2, int cj1, int cj2, amp::ampf< Precision > beta, ap::template_1d_array< amp::ampf< Precision > > &work)
int rowidxabsmax(const ap::template_2d_array< amp::ampf< Precision > > &x, int j1, int j2, int i)
void copyandtranspose(const ap::template_2d_array< amp::ampf< Precision > > &a, int is1, int is2, int js1, int js2, ap::template_2d_array< amp::ampf< Precision > > &b, int id1, int id2, int jd1, int jd2)
amp::ampf< Precision > vectornorm2(const ap::template_1d_array< amp::ampf< Precision > > &x, int i1, int i2)
amp::ampf< Precision > upperhessenberg1norm(const ap::template_2d_array< amp::ampf< Precision > > &a, int i1, int i2, int j1, int j2, ap::template_1d_array< amp::ampf< Precision > > &work)
void matrixvectormultiply(const ap::template_2d_array< amp::ampf< Precision > > &a, int i1, int i2, int j1, int j2, bool trans, const ap::template_1d_array< amp::ampf< Precision > > &x, int ix1, int ix2, amp::ampf< Precision > alpha, ap::template_1d_array< amp::ampf< Precision > > &y, int iy1, int iy2, amp::ampf< Precision > beta)
amp::ampf< Precision > pythag2(amp::ampf< Precision > x, amp::ampf< Precision > y)
void inplacetranspose(ap::template_2d_array< amp::ampf< Precision > > &a, int i1, int i2, int j1, int j2, ap::template_1d_array< amp::ampf< Precision > > &work)
void copymatrix(const ap::template_2d_array< amp::ampf< Precision > > &a, int is1, int is2, int js1, int js2, ap::template_2d_array< amp::ampf< Precision > > &b, int id1, int id2, int jd1, int jd2)
int vectoridxabsmax(const ap::template_1d_array< amp::ampf< Precision > > &x, int i1, int i2)
void rmatrixlq(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
void rmatrixlqunpackl(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &l)
void lqdecomposition(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
void rmatrixlqunpackq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qrows, ap::template_2d_array< amp::ampf< Precision > > &q)
void unpackqfromlq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qrows, ap::template_2d_array< amp::ampf< Precision > > &q)
void lqdecompositionunpacked(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &l, ap::template_2d_array< amp::ampf< Precision > > &q)
void rmatrixqr(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
void rmatrixqrunpackq(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
void qrdecomposition(ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_1d_array< amp::ampf< Precision > > &tau)
void unpackqfromqr(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, const ap::template_1d_array< amp::ampf< Precision > > &tau, int qcolumns, ap::template_2d_array< amp::ampf< Precision > > &q)
void rmatrixqrunpackr(const ap::template_2d_array< amp::ampf< Precision > > &a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &r)
void qrdecompositionunpacked(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, ap::template_2d_array< amp::ampf< Precision > > &q, ap::template_2d_array< amp::ampf< Precision > > &r)
void generatereflection(ap::template_1d_array< amp::ampf< Precision > > &x, int n, amp::ampf< Precision > &tau)
void applyreflectionfromtheright(ap::template_2d_array< amp::ampf< Precision > > &c, amp::ampf< Precision > tau, const ap::template_1d_array< amp::ampf< Precision > > &v, int m1, int m2, int n1, int n2, ap::template_1d_array< amp::ampf< Precision > > &work)
void applyreflectionfromtheleft(ap::template_2d_array< amp::ampf< Precision > > &c, amp::ampf< Precision > tau, const ap::template_1d_array< amp::ampf< Precision > > &v, int m1, int m2, int n1, int n2, ap::template_1d_array< amp::ampf< Precision > > &work)
void applyrotationsfromtheright(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
void generaterotation(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > &cs, amp::ampf< Precision > &sn, amp::ampf< Precision > &r)
void applyrotationsfromtheleft(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
bool rmatrixsvd(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)
bool svddecomposition(ap::template_2d_array< amp::ampf< Precision > > a, int m, int n, int uneeded, int vtneeded, int additionalmemory, ap::template_1d_array< amp::ampf< Precision > > &w, ap::template_2d_array< amp::ampf< Precision > > &u, ap::template_2d_array< amp::ampf< Precision > > &vt)