@@ -358,7 +358,7 @@ svd<Matrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
358
358
GESVD_REAL_STEP (dgesvd, DGESVD);
359
359
360
360
lwork = static_cast <F77_INT> (work[0 ]);
361
- work.reserve (lwork);
361
+ work.resize (lwork);
362
362
363
363
GESVD_REAL_STEP (dgesvd, DGESVD);
364
364
}
@@ -375,7 +375,7 @@ svd<FloatMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
375
375
GESVD_REAL_STEP (sgesvd, SGESVD);
376
376
377
377
lwork = static_cast <F77_INT> (work[0 ]);
378
- work.reserve (lwork);
378
+ work.resize (lwork);
379
379
380
380
GESVD_REAL_STEP (sgesvd, SGESVD);
381
381
}
@@ -394,7 +394,7 @@ svd<ComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
394
394
GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
395
395
396
396
lwork = static_cast <F77_INT> (work[0 ].real ());
397
- work.reserve (lwork);
397
+ work.resize (lwork);
398
398
399
399
GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
400
400
}
@@ -414,7 +414,7 @@ svd<FloatComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m,
414
414
GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
415
415
416
416
lwork = static_cast <F77_INT> (work[0 ].real ());
417
- work.reserve (lwork);
417
+ work.resize (lwork);
418
418
419
419
GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
420
420
}
@@ -450,7 +450,7 @@ svd<Matrix>::gesdd (char& jobz, F77_INT m, F77_INT n, double *tmp_data,
450
450
GESDD_REAL_STEP (dgesdd, DGESDD);
451
451
452
452
lwork = static_cast <F77_INT> (work[0 ]);
453
- work.reserve (lwork);
453
+ work.resize (lwork);
454
454
455
455
GESDD_REAL_STEP (dgesdd, DGESDD);
456
456
}
@@ -466,7 +466,7 @@ svd<FloatMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n, float *tmp_data,
466
466
GESDD_REAL_STEP (sgesdd, SGESDD);
467
467
468
468
lwork = static_cast <F77_INT> (work[0 ]);
469
- work.reserve (lwork);
469
+ work.resize (lwork);
470
470
471
471
GESDD_REAL_STEP (sgesdd, SGESDD);
472
472
}
@@ -495,7 +495,7 @@ svd<ComplexMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n,
495
495
GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
496
496
497
497
lwork = static_cast <F77_INT> (work[0 ].real ());
498
- work.reserve (lwork);
498
+ work.resize (lwork);
499
499
500
500
GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
501
501
}
@@ -524,7 +524,7 @@ svd<FloatComplexMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n,
524
524
GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
525
525
526
526
lwork = static_cast <F77_INT> (work[0 ].real ());
527
- work.reserve (lwork);
527
+ work.resize (lwork);
528
528
529
529
GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
530
530
}
@@ -581,7 +581,7 @@ svd<Matrix>::gejsv (char& joba, char& jobu, char& jobv,
581
581
F77_INT& info)
582
582
{
583
583
lwork = gejsv_lwork<Matrix>::optimal (joba, jobu, jobv, m, n);
584
- work.reserve (lwork);
584
+ work.resize (lwork);
585
585
586
586
GEJSV_REAL_STEP (dgejsv, DGEJSV);
587
587
}
@@ -598,7 +598,7 @@ svd<FloatMatrix>::gejsv (char& joba, char& jobu, char& jobv,
598
598
F77_INT& info)
599
599
{
600
600
lwork = gejsv_lwork<FloatMatrix>::optimal (joba, jobu, jobv, m, n);
601
- work.reserve (lwork);
601
+ work.resize (lwork);
602
602
603
603
GEJSV_REAL_STEP (sgejsv, SGEJSV);
604
604
}
@@ -616,18 +616,18 @@ svd<ComplexMatrix>::gejsv (char& joba, char& jobu, char& jobv,
616
616
{
617
617
F77_INT lrwork = -1 ; // work space size query
618
618
std::vector<double > rwork (1 );
619
- work.reserve (2 );
619
+ work.resize (2 );
620
620
621
621
GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
622
622
623
623
lwork = static_cast <F77_INT> (work[0 ].real ());
624
- work.reserve (lwork);
624
+ work.resize (lwork);
625
625
626
626
lrwork = static_cast <F77_INT> (rwork[0 ]);
627
- rwork.reserve (lrwork);
627
+ rwork.resize (lrwork);
628
628
629
629
F77_INT liwork = static_cast <F77_INT> (iwork[0 ]);
630
- iwork.reserve (liwork);
630
+ iwork.resize (liwork);
631
631
632
632
GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
633
633
}
@@ -645,18 +645,18 @@ svd<FloatComplexMatrix>::gejsv (char& joba, char& jobu, char& jobv,
645
645
{
646
646
F77_INT lrwork = -1 ; // work space size query
647
647
std::vector<float > rwork (1 );
648
- work.reserve (2 );
648
+ work.resize (2 );
649
649
650
650
GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
651
651
652
652
lwork = static_cast <F77_INT> (work[0 ].real ());
653
- work.reserve (lwork);
653
+ work.resize (lwork);
654
654
655
655
lrwork = static_cast <F77_INT> (rwork[0 ]);
656
- rwork.reserve (lrwork);
656
+ rwork.resize (lrwork);
657
657
658
658
F77_INT liwork = static_cast <F77_INT> (iwork[0 ]);
659
- iwork.reserve (liwork);
659
+ iwork.resize (liwork);
660
660
661
661
GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
662
662
}
@@ -758,9 +758,18 @@ svd<T>::svd (const T& a, svd::Type type, svd::Driver driver)
758
758
P *vt = m_right_sm.rwdata ();
759
759
760
760
// Query _GESVD for the correct dimension of WORK.
761
-
762
761
F77_INT lwork = -1 ;
763
-
762
+ // FIXME: A variable-sized scratchpad is required for Fortran SVD routines.
763
+ // Octave calls LAPACK routines with lwork of -1 initially which causes
764
+ // LAPACK to calculate the size of the required scratchpad and return that
765
+ // value in the first entry of the "work" array. The temporary buffer was
766
+ // grown to the desired size using std::vector::reserve(). However, this
767
+ // leads to undefined behavior outside the C++ specification. To work around
768
+ // this (7/4//2025) all of the calls were changed to std::vector::resize().
769
+ // This is correct, but unnecessarily initializes all values of the array.
770
+ // The code should be re-architected to avoid this, possibly by moving
771
+ // this scratchpad memory out of this function and in to the lower-level
772
+ // routines where the memory is used.
764
773
std::vector<P> work (1 );
765
774
766
775
const F77_INT f77_int_one = static_cast <F77_INT> (1 );
0 commit comments