diff --git a/cmake/lapack.cmake b/cmake/lapack.cmake index 04dd49cfbe..316aee4c81 100644 --- a/cmake/lapack.cmake +++ b/cmake/lapack.cmake @@ -11,8 +11,9 @@ set(SCLAUX la_constants.f90 sbdsdc.f sbdsqr.f sdisna.f slabad.f slacpy.f sladiv.f slae2.f slaebz.f - slaed0.f slaed1.f slaed2.f slaed4.f slaed5.f slaed6.f - slaed7.f slaed8.f slaed9.f slaeda.f slaev2.f slagtf.f + slaed0.f slaed0_br.f slaed1.f slaed2.f slaed4.f slaed4_br.f slaed5.f + slaed6.f slaed7.f slaed7_br.f slaed8.f slaed9.f slaed9_br.f + slaeda.f slaev2.f slagtf.f slagts.f slamrg.f slanst.f slapy2.f slapy3.f slarnv.f slarra.f slarrb.f slarrc.f slarrd.f slarre.f slarrf.f slarrj.f @@ -22,7 +23,7 @@ set(SCLAUX slasd7.f slasd8.f slasda.f slasdq.f slasdt.f slaset.f slasq1.f slasq2.f slasq3.f slasq4.f slasq5.f slasq6.f slasr.f slasrt.f slassq.f90 slasv2.f spttrf.f sstebz.f sstedc.f - ssteqr.f ssterf.f slaisnan.f sisnan.f slarmm.f + sstebr.f ssteqr.f ssterf.f slaisnan.f sisnan.f slarmm.f slartgp.f slartgs.f ../INSTALL/sroundup_lwork.f ../INSTALL/second_${TIMER}.f) @@ -31,8 +32,9 @@ set(DZLAUX dbdsdc.f dbdsvdx.f dstevx.f dstein.f dbdsqr.f ddisna.f dlabad.f dlacpy.f dladiv.f dlae2.f dlaebz.f - dlaed0.f dlaed1.f dlaed2.f dlaed4.f dlaed5.f dlaed6.f - dlaed7.f dlaed8.f dlaed9.f dlaeda.f dlaev2.f dlagtf.f + dlaed0.f dlaed0_br.f dlaed1.f dlaed2.f dlaed4.f dlaed4_br.f dlaed5.f + dlaed6.f dlaed7.f dlaed7_br.f dlaed8.f dlaed9.f dlaed9_br.f + dlaeda.f dlaev2.f dlagtf.f dlagts.f dlamrg.f dlanst.f dlapy2.f dlapy3.f dlarnv.f dlarra.f dlarrb.f dlarrc.f dlarrd.f dlarre.f dlarrf.f dlarrj.f @@ -42,7 +44,7 @@ set(DZLAUX dlasd7.f dlasd8.f dlasda.f dlasdq.f dlasdt.f dlasq1.f dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f dlasr.f dlasrt.f dlassq.f90 dlasv2.f dpttrf.f dstebz.f dstedc.f - dsteqr.f dsterf.f dlaisnan.f disnan.f dlarmm.f + dstebr.f dsteqr.f dsterf.f dlaisnan.f disnan.f dlarmm.f dlartgp.f dlartgs.f ../INSTALL/droundup_lwork.f ../INSTALL/dlamch.f ../INSTALL/dsecnd_${TIMER}.f) diff --git a/exports/gensymbol b/exports/gensymbol index 6284e6b6ae..9cb966007f 100755 --- a/exports/gensymbol +++ b/exports/gensymbol @@ -280,8 +280,8 @@ lapackobjs2sc=" slasd0 slasd1 slasd2 slasd3 slasd4 slasd5 slasd6 slasd7 slasd8 slasda slasdq slasdt slaset slasq1 slasq2 slasq3 slasq4 slasq5 slasq6 - slasr slasrt slassq slasv2 spttrf sstebz sstedc - ssteqr ssterf slaisnan sisnan + slasr slasrt slassq slasv2 spttrf sstebz sstebr + sstedc ssteqr ssterf slaisnan sisnan slartgp slartgs " @@ -300,8 +300,8 @@ lapackobjs2dz=" dlasd0 dlasd1 dlasd2 dlasd3 dlasd4 dlasd5 dlasd6 dlasd7 dlasd8 dlasda dlasdq dlasdt dlaset dlasq1 dlasq2 dlasq3 dlasq4 dlasq5 dlasq6 - dlasr dlasrt dlassq dlasv2 dpttrf dstebz dstedc - dsteqr dsterf dlaisnan disnan + dlasr dlasrt dlassq dlasv2 dpttrf dstebz dstebr + dstedc dsteqr dsterf dlaisnan disnan dlartgp dlartgs dlamch dlamc3 " @@ -2117,6 +2117,8 @@ lapackeobjsd=" LAPACKE_dsptrs_work LAPACKE_dstebz LAPACKE_dstebz_work + LAPACKE_dstebr + LAPACKE_dstebr_work LAPACKE_dstedc LAPACKE_dstedc_work LAPACKE_dstegr @@ -2736,6 +2738,8 @@ lapackeobjss=" LAPACKE_ssptrs_work LAPACKE_sstebz LAPACKE_sstebz_work + LAPACKE_sstebr + LAPACKE_sstebr_work LAPACKE_sstedc LAPACKE_sstedc_work LAPACKE_sstegr diff --git a/lapack-netlib/LAPACKE/include/lapack.h b/lapack-netlib/LAPACKE/include/lapack.h index 0ed9ad01a9..78f7b39147 100644 --- a/lapack-netlib/LAPACKE/include/lapack.h +++ b/lapack-netlib/LAPACKE/include/lapack.h @@ -16368,6 +16368,24 @@ void LAPACK_sstebz_base( #define LAPACK_sstebz(...) LAPACK_sstebz_base(__VA_ARGS__) #endif +#define LAPACK_dstebr LAPACK_GLOBAL(dstebr,DSTEBR) +void LAPACK_dstebr( + lapack_int const* n, + double* D, + double* E, + double* work, lapack_int const* lwork, + lapack_int* iwork, lapack_int const* liwork, + lapack_int* info ); + +#define LAPACK_sstebr LAPACK_GLOBAL(sstebr,SSTEBR) +void LAPACK_sstebr( + lapack_int const* n, + float* D, + float* E, + float* work, lapack_int const* lwork, + lapack_int* iwork, lapack_int const* liwork, + lapack_int* info ); + #define LAPACK_cstedc_base LAPACK_GLOBAL(cstedc,CSTEDC) void LAPACK_cstedc_base( char const* compz, diff --git a/lapack-netlib/LAPACKE/include/lapacke.h b/lapack-netlib/LAPACKE/include/lapacke.h index 377e2a6bbc..a14f347bc0 100644 --- a/lapack-netlib/LAPACKE/include/lapacke.h +++ b/lapack-netlib/LAPACKE/include/lapacke.h @@ -3574,6 +3574,9 @@ lapack_int LAPACKE_dstebz( char range, char order, lapack_int n, double vl, lapack_int* m, lapack_int* nsplit, double* w, lapack_int* iblock, lapack_int* isplit ); +lapack_int LAPACKE_sstebr( lapack_int n, float* d, float* e ); +lapack_int LAPACKE_dstebr( lapack_int n, double* d, double* e ); + lapack_int LAPACKE_sstedc( int matrix_layout, char compz, lapack_int n, float* d, float* e, float* z, lapack_int ldz ); lapack_int LAPACKE_dstedc( int matrix_layout, char compz, lapack_int n, @@ -9197,6 +9200,13 @@ lapack_int LAPACKE_dstebz_work( char range, char order, lapack_int n, double vl, lapack_int* iblock, lapack_int* isplit, double* work, lapack_int* iwork ); +lapack_int LAPACKE_sstebr_work( lapack_int n, float* d, float* e, float* work, + lapack_int lwork, lapack_int* iwork, + lapack_int liwork ); +lapack_int LAPACKE_dstebr_work( lapack_int n, double* d, double* e, + double* work, lapack_int lwork, + lapack_int* iwork, lapack_int liwork ); + lapack_int LAPACKE_sstedc_work( int matrix_layout, char compz, lapack_int n, float* d, float* e, float* z, lapack_int ldz, float* work, lapack_int lwork, diff --git a/lapack-netlib/LAPACKE/src/CMakeLists.txt b/lapack-netlib/LAPACKE/src/CMakeLists.txt index 89890f9239..9fb17efe2e 100644 --- a/lapack-netlib/LAPACKE/src/CMakeLists.txt +++ b/lapack-netlib/LAPACKE/src/CMakeLists.txt @@ -1011,6 +1011,8 @@ lapacke_dsptrs.c lapacke_dsptrs_work.c lapacke_dstebz.c lapacke_dstebz_work.c +lapacke_dstebr.c +lapacke_dstebr_work.c lapacke_dstedc.c lapacke_dstedc_work.c lapacke_dstegr.c @@ -1584,6 +1586,8 @@ lapacke_ssptrs.c lapacke_ssptrs_work.c lapacke_sstebz.c lapacke_sstebz_work.c +lapacke_sstebr.c +lapacke_sstebr_work.c lapacke_sstedc.c lapacke_sstedc_work.c lapacke_sstegr.c diff --git a/lapack-netlib/LAPACKE/src/Makefile b/lapack-netlib/LAPACKE/src/Makefile index 969288f424..128a060903 100644 --- a/lapack-netlib/LAPACKE/src/Makefile +++ b/lapack-netlib/LAPACKE/src/Makefile @@ -1075,6 +1075,8 @@ lapacke_dsptrs.o \ lapacke_dsptrs_work.o \ lapacke_dstebz.o \ lapacke_dstebz_work.o \ +lapacke_dstebr.o \ +lapacke_dstebr_work.o \ lapacke_dstedc.o \ lapacke_dstedc_work.o \ lapacke_dstegr.o \ @@ -1653,6 +1655,8 @@ lapacke_ssptrs.o \ lapacke_ssptrs_work.o \ lapacke_sstebz.o \ lapacke_sstebz_work.o \ +lapacke_sstebr.o \ +lapacke_sstebr_work.o \ lapacke_sstedc.o \ lapacke_sstedc_work.o \ lapacke_sstegr.o \ diff --git a/lapack-netlib/LAPACKE/src/lapacke_dstebr.c b/lapack-netlib/LAPACKE/src/lapacke_dstebr.c new file mode 100644 index 0000000000..0864526855 --- /dev/null +++ b/lapack-netlib/LAPACKE/src/lapacke_dstebr.c @@ -0,0 +1,87 @@ +/***************************************************************************** + Copyright (c) 2014, Intel Corp. + Copyright (c) 2026, Ruiyi Zhan and Shaoshuai Zhang. + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF + THE POSSIBILITY OF SUCH DAMAGE. +***************************************************************************** +* Contents: Native high-level C interface to LAPACK function dstebr +* Author: Ruiyi Zhan and Shaoshuai Zhang +*****************************************************************************/ + +#include "lapacke_utils.h" + +lapack_int LAPACKE_dstebr( lapack_int n, double* d, double* e ) +{ + lapack_int info = 0; + lapack_int liwork = -1; + lapack_int lwork = -1; + lapack_int* iwork = NULL; + double* work = NULL; + lapack_int iwork_query; + double work_query; +#ifndef LAPACK_DISABLE_NAN_CHECK + if( LAPACKE_get_nancheck() ) { + /* Optionally check input matrices for NaNs */ + if( LAPACKE_d_nancheck( n, d, 1 ) ) { + return -2; + } + if( LAPACKE_d_nancheck( n-1, e, 1 ) ) { + return -3; + } + } +#endif + /* Query optimal working array(s) size */ + info = LAPACKE_dstebr_work( n, d, e, &work_query, lwork, + &iwork_query, liwork ); + if( info != 0 ) { + goto exit_level_0; + } + liwork = iwork_query; + lwork = (lapack_int)work_query; + /* Allocate memory for work arrays */ + iwork = (lapack_int*)LAPACKE_malloc( sizeof(lapack_int) * liwork ); + if( iwork == NULL ) { + info = LAPACK_WORK_MEMORY_ERROR; + goto exit_level_0; + } + work = (double*)LAPACKE_malloc( sizeof(double) * lwork ); + if( work == NULL ) { + info = LAPACK_WORK_MEMORY_ERROR; + goto exit_level_1; + } + /* Call middle-level interface */ + info = LAPACKE_dstebr_work( n, d, e, work, lwork, iwork, + liwork ); + /* Release memory and exit */ + LAPACKE_free( work ); +exit_level_1: + LAPACKE_free( iwork ); +exit_level_0: + if( info == LAPACK_WORK_MEMORY_ERROR ) { + LAPACKE_xerbla( "LAPACKE_dstebr", info ); + } + return info; +} diff --git a/lapack-netlib/LAPACKE/src/lapacke_dstebr_work.c b/lapack-netlib/LAPACKE/src/lapacke_dstebr_work.c new file mode 100644 index 0000000000..d71c2d5317 --- /dev/null +++ b/lapack-netlib/LAPACKE/src/lapacke_dstebr_work.c @@ -0,0 +1,44 @@ +/***************************************************************************** + Copyright (c) 2014, Intel Corp. + Copyright (c) 2026, Ruiyi Zhan and Shaoshuai Zhang. + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF + THE POSSIBILITY OF SUCH DAMAGE. +***************************************************************************** +* Contents: Native middle-level C interface to LAPACK function dstebr +* Author: Ruiyi Zhan and Shaoshuai Zhang +*****************************************************************************/ + +#include "lapacke_utils.h" + +lapack_int LAPACKE_dstebr_work( lapack_int n, double* d, double* e, + double* work, lapack_int lwork, + lapack_int* iwork, lapack_int liwork ) +{ + lapack_int info = 0; + /* Call LAPACK function and adjust info */ + LAPACK_dstebr( &n, d, e, work, &lwork, iwork, &liwork, &info ); + return info; +} diff --git a/lapack-netlib/LAPACKE/src/lapacke_sstebr.c b/lapack-netlib/LAPACKE/src/lapacke_sstebr.c new file mode 100644 index 0000000000..67cfce2ee9 --- /dev/null +++ b/lapack-netlib/LAPACKE/src/lapacke_sstebr.c @@ -0,0 +1,87 @@ +/***************************************************************************** + Copyright (c) 2014, Intel Corp. + Copyright (c) 2026, Ruiyi Zhan and Shaoshuai Zhang. + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF + THE POSSIBILITY OF SUCH DAMAGE. +***************************************************************************** +* Contents: Native high-level C interface to LAPACK function sstebr +* Author: Ruiyi Zhan and Shaoshuai Zhang +*****************************************************************************/ + +#include "lapacke_utils.h" + +lapack_int LAPACKE_sstebr( lapack_int n, float* d, float* e ) +{ + lapack_int info = 0; + lapack_int liwork = -1; + lapack_int lwork = -1; + lapack_int* iwork = NULL; + float* work = NULL; + lapack_int iwork_query; + float work_query; +#ifndef LAPACK_DISABLE_NAN_CHECK + if( LAPACKE_get_nancheck() ) { + /* Optionally check input matrices for NaNs */ + if( LAPACKE_s_nancheck( n, d, 1 ) ) { + return -2; + } + if( LAPACKE_s_nancheck( n-1, e, 1 ) ) { + return -3; + } + } +#endif + /* Query optimal working array(s) size */ + info = LAPACKE_sstebr_work( n, d, e, &work_query, lwork, + &iwork_query, liwork ); + if( info != 0 ) { + goto exit_level_0; + } + liwork = iwork_query; + lwork = (lapack_int)work_query; + /* Allocate memory for work arrays */ + iwork = (lapack_int*)LAPACKE_malloc( sizeof(lapack_int) * liwork ); + if( iwork == NULL ) { + info = LAPACK_WORK_MEMORY_ERROR; + goto exit_level_0; + } + work = (float*)LAPACKE_malloc( sizeof(float) * lwork ); + if( work == NULL ) { + info = LAPACK_WORK_MEMORY_ERROR; + goto exit_level_1; + } + /* Call middle-level interface */ + info = LAPACKE_sstebr_work( n, d, e, work, lwork, iwork, + liwork ); + /* Release memory and exit */ + LAPACKE_free( work ); +exit_level_1: + LAPACKE_free( iwork ); +exit_level_0: + if( info == LAPACK_WORK_MEMORY_ERROR ) { + LAPACKE_xerbla( "LAPACKE_sstebr", info ); + } + return info; +} diff --git a/lapack-netlib/LAPACKE/src/lapacke_sstebr_work.c b/lapack-netlib/LAPACKE/src/lapacke_sstebr_work.c new file mode 100644 index 0000000000..738b08935f --- /dev/null +++ b/lapack-netlib/LAPACKE/src/lapacke_sstebr_work.c @@ -0,0 +1,44 @@ +/***************************************************************************** + Copyright (c) 2014, Intel Corp. + Copyright (c) 2026, Ruiyi Zhan and Shaoshuai Zhang. + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF + THE POSSIBILITY OF SUCH DAMAGE. +***************************************************************************** +* Contents: Native middle-level C interface to LAPACK function sstebr +* Author: Ruiyi Zhan and Shaoshuai Zhang +*****************************************************************************/ + +#include "lapacke_utils.h" + +lapack_int LAPACKE_sstebr_work( lapack_int n, float* d, float* e, + float* work, lapack_int lwork, + lapack_int* iwork, lapack_int liwork ) +{ + lapack_int info = 0; + /* Call LAPACK function and adjust info */ + LAPACK_sstebr( &n, d, e, work, &lwork, iwork, &liwork, &info ); + return info; +} diff --git a/lapack-netlib/SRC/CMakeLists.txt b/lapack-netlib/SRC/CMakeLists.txt index 90c5e61a98..be8007ad9c 100644 --- a/lapack-netlib/SRC/CMakeLists.txt +++ b/lapack-netlib/SRC/CMakeLists.txt @@ -44,8 +44,10 @@ set(SCLAUX la_constants.f90 sbdsvdx.f sbdsdc.f sbdsqr.f sdisna.f slabad.f slacpy.f sladiv.f slae2.f slaebz.f - slaed0.f slaed1.f slaed2.f slaed3.f slaed4.f slaed5.f slaed6.f - slaed7.f slaed8.f slaed9.f slaeda.f slaev2.f slagtf.f + slaed0.f slaed0_br.f slaed1.f slaed2.f slaed3.f slaed4.f + slaed4_br.f slaed5.f slaed6.f + slaed7.f slaed7_br.f slaed8.f slaed9.f slaed9_br.f slaeda.f + slaev2.f slagtf.f slagts.f slamrg.f slanst.f slapy2.f slapy3.f slarnv.f slarra.f slarrb.f slarrc.f slarrd.f slarre.f slarrf.f slarrj.f @@ -55,7 +57,7 @@ set(SCLAUX slasd7.f slasd8.f slasda.f slasdq.f slasdt.f slaset.f slasq1.f slasq2.f slasq3.f slasq4.f slasq5.f slasq6.f slasr.f slasrt.f slassq.f90 slasv2.f spttrf.f sstebz.f sstedc.f - sstein.f ssteqr.f ssterf.f sstevx.f + sstein.f sstebr.f ssteqr.f ssterf.f sstevx.f slartgp.f slartgs.f ../INSTALL/sroundup_lwork.f ${SECOND_SRC}) @@ -65,8 +67,10 @@ set(DZLAUX dbdsqr.f ddisna.f disnan.f dlabad.f dlacpy.f dladiv.f dlae2.f dlaebz.f - dlaed0.f dlaed1.f dlaed2.f dlaed3.f dlaed4.f dlaed5.f dlaed6.f - dlaed7.f dlaed8.f dlaed9.f dlaeda.f dlaev2.f dlagtf.f + dlaed0.f dlaed0_br.f dlaed1.f dlaed2.f dlaed3.f dlaed4.f + dlaed4_br.f dlaed5.f dlaed6.f + dlaed7.f dlaed7_br.f dlaed8.f dlaed9.f dlaed9_br.f dlaeda.f + dlaev2.f dlagtf.f dlagts.f dlamrg.f dlanst.f dlapy2.f dlapy3.f dlarnv.f dlarra.f dlarrb.f dlarrc.f dlarrd.f dlarre.f dlarrf.f dlarrj.f @@ -77,7 +81,7 @@ set(DZLAUX dlaset.f dlasq1.f dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f dlasr.f dlasrt.f dlassq.f90 dlasv2.f dlaisnan.f dpttrf.f - dstebz.f dstedc.f dstein.f dsteqr.f dsterf.f dstevx.f + dstebz.f dstedc.f dstebr.f dstein.f dsteqr.f dsterf.f dstevx.f dlartgp.f dlartgs.f ../INSTALL/droundup_lwork.f ../INSTALL/dlamch.f ${DSECOND_SRC}) diff --git a/lapack-netlib/SRC/Makefile b/lapack-netlib/SRC/Makefile index ebf3431a92..74cff4ad88 100644 --- a/lapack-netlib/SRC/Makefile +++ b/lapack-netlib/SRC/Makefile @@ -103,6 +103,9 @@ SCLAUX_O = la_constants.o \ ssteqr.o ssterf.o slaisnan.o sisnan.o \ slartgp.o slartgs.o scombssq.o ../INSTALL/sroundup_lwork.o \ ../INSTALL/second_$(TIMER).o +ifneq ($(C_LAPACK), 1) +SCLAUX_O += slaed0_br.o slaed4_br.o slaed7_br.o slaed9_br.o sstebr.o +endif endif ifneq "$(or $(BUILD_DOUBLE),$(BUILD_COMPLEX16))" "" @@ -125,6 +128,9 @@ DZLAUX_O = la_constants.o\ dsteqr.o dsterf.o dlaisnan.o disnan.o \ dlartgp.o dlartgs.o ../INSTALL/droundup_lwork.o \ ../INSTALL/dlamch.o ../INSTALL/dsecnd_$(TIMER).o +ifneq ($(C_LAPACK), 1) +DZLAUX_O += dlaed0_br.o dlaed4_br.o dlaed7_br.o dlaed9_br.o dstebr.o +endif endif #ifeq ($(BUILD_SINGLE),1) diff --git a/lapack-netlib/SRC/dlaed0_br.f b/lapack-netlib/SRC/dlaed0_br.f new file mode 100644 index 0000000000..98b8de0592 --- /dev/null +++ b/lapack-netlib/SRC/dlaed0_br.f @@ -0,0 +1,233 @@ +*> \brief \b DLAED0_BR used by DSTEBR. +* +* DLAED0_BR computes eigenvalues of a symmetric tridiagonal matrix with +* the boundary-row divide-and-conquer method used by DSTEBR. Unlike +* DLAED0(ICOMPQ=0), it keeps only first/last boundary rows for each +* active block across merge levels instead of saving complete +* eigenvector blocks. +* +* Serial Reference-LAPACK form: no OpenMP, environment-variable +* scheduling, or thread-scaled workspace. +* + SUBROUTINE DLAED0_BR( N, D, E, WORK, IWORK, INFO ) +* +* .. Scalar Arguments .. + INTEGER INFO, N +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION D( * ), E( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + INTEGER WRKSTR + PARAMETER ( WRKSTR = 14 ) +* .. +* .. Local Scalars .. + LOGICAL WANTQ + INTEGER I, IBHI, IBLO, IIWRK, INDXQ, INFOM, ISCR, + $ IWBASE, J, K, KDEFL, MATSIZ, MERGE, MSD2, + $ NMERGE, QBASE, SMLSIZ, SMM1, SPM1, SUBMAT, + $ SUBPBS, WBASE, WRKBASE +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DLAED7_BR, DSTEQR_BOUNDARY_BR, XERBLA +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( N.LT.0 ) THEN + INFO = -1 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAED0_BR', -INFO ) + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN +* + SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 ) +* +* Determine the size and placement of leaf submatrices. +* + IWORK( 1 ) = N + SUBPBS = 1 + 10 CONTINUE + IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN + DO 20 J = SUBPBS, 1, -1 + IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 + IWORK( 2*J-1 ) = IWORK( J ) / 2 + 20 CONTINUE + SUBPBS = 2*SUBPBS + GO TO 10 + END IF + DO 30 J = 2, SUBPBS + IWORK( J ) = IWORK( J ) + IWORK( J-1 ) + 30 CONTINUE +* +* Apply rank-1 cuts. +* + SPM1 = SUBPBS - 1 + DO 40 I = 1, SPM1 + SUBMAT = IWORK( I ) + 1 + SMM1 = SUBMAT - 1 + D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) + D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) + 40 CONTINUE +* +* Workspace layout: +* WORK(IBLO:IBLO+N-1) first boundary rows +* WORK(IBHI:IBHI+N-1) last boundary rows +* WORK(ISCR:) leaf scratch and per-merge scratch +* +* Integer workspace layout: +* IWORK(1:N) split endpoints +* IWORK(INDXQ:INDXQ+N-1) sorting permutations +* IWORK(IIWRK:) DLAED7_BR scratch +* + IBLO = 1 + IBHI = IBLO + N + ISCR = IBHI + N + INDXQ = N + 1 + IIWRK = INDXQ + N +* +* Solve leaf eigenproblems and initialize boundary row invariants. +* + DO 70 I = 1, SUBPBS + QBASE = ISCR + WBASE = QBASE + SMLSIZ*SMLSIZ + IF( I.EQ.1 ) THEN + SUBMAT = 1 + MATSIZ = IWORK( 1 ) + ELSE + SUBMAT = IWORK( I-1 ) + 1 + MATSIZ = IWORK( I ) - IWORK( I-1 ) + END IF +* + CALL DSTEQR_BOUNDARY_BR( MATSIZ, D( SUBMAT ), E( SUBMAT ), + $ WORK( IBLO+SUBMAT-1 ), + $ WORK( IBHI+SUBMAT-1 ), + $ WORK( QBASE ), SMLSIZ, WORK( WBASE ), + $ INFO ) + IF( INFO.NE.0 ) + $ GO TO 130 + K = 1 + DO 65 J = SUBMAT, IWORK( I ) + IWORK( INDXQ+J-1 ) = K + K = K + 1 + 65 CONTINUE + 70 CONTINUE +* +* Merge adjacent eigensystems. Boundary rows are updated in place +* for the merged physical D order. +* + 80 CONTINUE + IF( SUBPBS.GT.1 ) THEN + WANTQ = SUBPBS.GT.2 + NMERGE = SUBPBS / 2 + DO 90 MERGE = 1, NMERGE + CALL DLAED0_MERGE_PAIR( MERGE, NMERGE, IWORK, WANTQ, D, E, + $ WORK, IBLO, IBHI, ISCR, IIWRK, INDXQ, SUBMAT, MATSIZ, + $ MSD2, WRKBASE, IWBASE, KDEFL, INFOM, WRKSTR ) + IF( INFOM.NE.0 ) THEN + INFO = INFOM + GO TO 130 + END IF + 90 CONTINUE +* + DO 95 MERGE = 1, NMERGE + IWORK( MERGE ) = IWORK( 2*MERGE ) + 95 CONTINUE + SUBPBS = SUBPBS / 2 + GO TO 80 + END IF +* +* Reorder eigenvalues from the final physical merge order. +* + DO 100 I = 1, N + J = IWORK( INDXQ+I-1 ) + WORK( ISCR+I-1 ) = D( J ) + 100 CONTINUE + CALL DCOPY( N, WORK( ISCR ), 1, D, 1 ) + GO TO 140 +* + 130 CONTINUE + INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 +* + 140 CONTINUE + RETURN +* +* End of DLAED0_BR +* + END +* +*> \brief \b DLAED0_MERGE_PAIR executes one merge in the DC tree. +* + SUBROUTINE DLAED0_MERGE_PAIR( MERGE, NMERGE, IWORK, WANTQ, D, E, + $ WORK, IBLO, IBHI, ISCR, IIWRK, + $ INDXQ, SUBMAT, MATSIZ, MSD2, + $ WRKBASE, IWBASE, KDEFL, INFO, + $ WRKSTR ) +* + INTEGER IBHI, IBLO, IIWRK, INDXQ, INFO, ISCR, IWBASE, + $ KDEFL, MATSIZ, MERGE, MSD2, NMERGE, RIGHT, + $ SUBMAT, WRKBASE, WRKSTR + LOGICAL WANTQ + INTEGER IWORK( * ) + DOUBLE PRECISION D( * ), E( * ), WORK( * ) + EXTERNAL DLAED7_BR +* + INTEGER LEFT +* + LEFT = 2*MERGE - 1 + RIGHT = LEFT + 1 + IF( MERGE.EQ.1 ) THEN + SUBMAT = 1 + MATSIZ = IWORK( RIGHT ) + MSD2 = IWORK( 1 ) + ELSE + SUBMAT = IWORK( LEFT-1 ) + 1 + MATSIZ = IWORK( RIGHT ) - IWORK( LEFT-1 ) + MSD2 = IWORK( LEFT ) - IWORK( LEFT-1 ) + END IF +* + WRKBASE = ISCR + WRKSTR*( SUBMAT-1 ) + IWBASE = IIWRK + 5*( SUBMAT-1 ) +* + CALL DLAED7_BR( WANTQ, MATSIZ, MSD2, D( SUBMAT ), + $ IWORK( INDXQ+SUBMAT-1 ), E( SUBMAT+MSD2-1 ), + $ WORK( IBLO+SUBMAT-1 ), WORK( IBHI+SUBMAT-1 ), + $ WORK( WRKBASE ), IWORK( IWBASE ), KDEFL, INFO ) + RETURN + END +* +*> \brief \b DSTEQR_BOUNDARY_BR solves one leaf and keeps boundary rows. +* + SUBROUTINE DSTEQR_BOUNDARY_BR( N, D, E, BLO, BHI, Q, LDQ, WORK, + $ INFO ) +* + INTEGER INFO, J, LDQ, N + DOUBLE PRECISION BLO( * ), BHI( * ), D( * ), E( * ), Q( * ), + $ WORK( * ) + EXTERNAL DSTEQR +* + CALL DSTEQR( 'I', N, D, E, Q, LDQ, WORK, INFO ) + IF( INFO.NE.0 ) + $ RETURN +* + DO 10 J = 1, N + BLO( J ) = Q( 1 + ( J - 1 )*LDQ ) + BHI( J ) = Q( N + ( J - 1 )*LDQ ) + 10 CONTINUE + RETURN + END diff --git a/lapack-netlib/SRC/dlaed4_br.f b/lapack-netlib/SRC/dlaed4_br.f new file mode 100644 index 0000000000..fec1e05a41 --- /dev/null +++ b/lapack-netlib/SRC/dlaed4_br.f @@ -0,0 +1,109 @@ +*> \brief \b DLAED4_BR used by DSTEBR. +* +* DLAED4_BR computes one secular root with DLAED4 and exposes a compact +* representation of the returned high-relative-accuracy DELTA vector. +* For N > 2, DLAED4's DELTA entries are all relative to either D(I) or +* D(I+1), except for the final root where the origin is D(N). ORG records +* that origin index and TAU satisfies +* +* lambda = D(ORG) + TAU, +* DELTA(j) = D(j) - D(ORG) - TAU. +* +* NEAR1/NEAR2 cache the two pole-adjacent DELTA entries that are most +* sensitive to cancellation when reconstructed from D and TAU. +* + SUBROUTINE DLAED4_BR( N, I, D, Z, DELTA, RHO, DLAM, ORG, + $ TAU, NEAR1, NEAR2, INFO ) +* +* .. Scalar Arguments .. + INTEGER I, INFO, N, ORG + DOUBLE PRECISION DLAM, NEAR1, NEAR2, RHO, TAU +* .. +* .. Array Arguments .. + DOUBLE PRECISION D( * ), DELTA( * ), Z( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION HALF, ZERO + PARAMETER ( HALF = 0.5D0, ZERO = 0.0D0 ) +* .. +* .. Local Scalars .. + DOUBLE PRECISION MID +* .. +* .. External Subroutines .. + EXTERNAL DLAED4, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* + INFO = 0 + ORG = 1 + TAU = ZERO + NEAR1 = ZERO + NEAR2 = ZERO +* + IF( N.LT.1 ) THEN + INFO = -1 + ELSE IF( I.LT.1 .OR. I.GT.N ) THEN + INFO = -2 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAED4_BR', -INFO ) + RETURN + END IF +* + CALL DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) + IF( INFO.NE.0 ) + $ RETURN +* + IF( N.EQ.1 ) THEN + ORG = 1 + TAU = DLAM - D( 1 ) + NEAR1 = DELTA( 1 ) + NEAR2 = ZERO + RETURN + END IF +* + IF( N.EQ.2 ) THEN + ORG = 1 + TAU = DLAM - D( 1 ) + NEAR1 = DELTA( 1 ) + NEAR2 = DELTA( 2 ) + RETURN + END IF +* + IF( I.EQ.N ) THEN + ORG = N + ELSE + MID = HALF*( D( I )+D( I+1 ) ) + IF( DLAM.LT.MID ) THEN + ORG = I + ELSE + ORG = I + 1 + END IF + END IF +* +* Avoid reconstructing TAU from DLAM-D(ORG). DLAED4 already +* returned DELTA(ORG) in its shifted coordinate, so using +* TAU=-DELTA(ORG) preserves that final rounded denominator and +* avoids an extra add/subtract rounding through DLAM. +* + TAU = -DELTA( ORG ) +* + IF( I.EQ.N ) THEN + NEAR1 = DELTA( N-1 ) + NEAR2 = DELTA( N ) + ELSE + NEAR1 = DELTA( I ) + NEAR2 = DELTA( I+1 ) + END IF +* + RETURN +* +* End of DLAED4_BR +* + END diff --git a/lapack-netlib/SRC/dlaed7_br.f b/lapack-netlib/SRC/dlaed7_br.f new file mode 100644 index 0000000000..0b47071fc7 --- /dev/null +++ b/lapack-netlib/SRC/dlaed7_br.f @@ -0,0 +1,376 @@ +*> \brief \b DLAED7_BR used by DSTEBR. +* +* DLAED7_BR performs one divide-and-conquer merge while propagating +* only the first and last rows of the current local eigenvector block. +* It is called by DLAED0_BR and is not itself a public LAPACK API. +* +* Serial Reference-LAPACK form: no OpenMP or thread-scaled workspace. +* + SUBROUTINE DLAED7_BR( WANTQ, N, CUTPNT, D, INDXQ, RHO, BLO, BHI, + $ WORK, IWORK, KDEFL, INFO ) +* +* .. Scalar Arguments .. + LOGICAL WANTQ + INTEGER CUTPNT, INFO, KDEFL, N + DOUBLE PRECISION RHO +* .. +* .. Array Arguments .. + INTEGER INDXQ( * ), IWORK( * ) + DOUBLE PRECISION BLO( * ), BHI( * ), D( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) +* .. +* .. Local Scalars .. + INTEGER GIVPTR, IDLMDA, IGIVCL, IGIVNM, I, IINDX, + $ IINDXP, INEAR, IPERM, IQ, IQ2, ISCR, ITAU, IW, + $ IWSGN, IZ, J, K, N1, N2, NBR +* .. +* .. External Subroutines .. + EXTERNAL DLAED8_BR, DLAED9_BR, DLAMRG, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* + INFO = 0 + KDEFL = 0 + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN + INFO = -2 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAED7_BR', -INFO ) + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN +* +* Workspace layout (14*N doubles): +* +* Z, DLAMBDA, W : N each +* QBR, Q2BR, GIVNUM/reused WPROD : 2*N each +* DELTA, WSGN, TAU, NEAR : N, N, N, 2*N +* + IZ = 1 + IDLMDA = IZ + N + IW = IDLMDA + N + IQ = IW + N + IQ2 = IQ + 2*N + IGIVNM = IQ2 + 2*N + ISCR = IGIVNM + 2*N + IWSGN = ISCR + N + ITAU = IWSGN + N + INEAR = ITAU + N +* +* Integer workspace layout. +* + IINDXP = 1 + IINDX = IINDXP + N + IPERM = IINDX + N + IGIVCL = IPERM + N +* + IF( WANTQ ) THEN + NBR = 2 + ELSE + NBR = 0 + END IF +* +* Z is the rank-one update vector. QBR has two selected rows of +* blockdiag(Q_left,Q_right): parent first boundary and parent last +* boundary, represented in the child eigenbases. The root merge has +* no parent, so it only needs Z for the secular eigenvalues. +* + IF( WANTQ ) THEN + DO 20 J = 1, CUTPNT + WORK( IZ+J-1 ) = BHI( J ) + WORK( IQ+2*( J-1 ) ) = BLO( J ) + WORK( IQ+2*( J-1 )+1 ) = ZERO + 20 CONTINUE + DO 30 J = CUTPNT + 1, N + WORK( IZ+J-1 ) = BLO( J ) + WORK( IQ+2*( J-1 ) ) = ZERO + WORK( IQ+2*( J-1 )+1 ) = BHI( J ) + 30 CONTINUE + ELSE + DO 40 J = 1, CUTPNT + WORK( IZ+J-1 ) = BHI( J ) + 40 CONTINUE + DO 50 J = CUTPNT + 1, N + WORK( IZ+J-1 ) = BLO( J ) + 50 CONTINUE + END IF +* +* Sort and deflate, applying the same Givens/permutation actions to +* the selected boundary rows when a parent merge will need them. +* +* + CALL DLAED8_BR( NBR, K, N, D, WORK( IQ ), 2, INDXQ, RHO, + $ CUTPNT, WORK( IZ ), WORK( IDLMDA ), + $ WORK( IQ2 ), 2, WORK( IW ), IWORK( IPERM ), + $ GIVPTR, IWORK( IGIVCL ), WORK( IGIVNM ), + $ IWORK( IINDXP ), IWORK( IINDX ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 70 + KDEFL = K +* +* Solve the secular equation and directly update the selected rows. +* WORK(ISCR) holds one DLAED4 delta vector at a time; no K-by-K +* delta or eigenvector block is materialized. +* + IF( K.NE.0 ) THEN + CALL DLAED9_BR( K, NBR, D, WORK( IQ2 ), 2, WORK( IQ ), 2, + $ RHO, WORK( IDLMDA ), WORK( IW ), + $ WORK( ISCR ), WORK( IGIVNM ), + $ WORK( IWSGN ), WORK( ITAU ), + $ WORK( INEAR ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 70 +* +* Prepare INDXQ sorting permutation. +* + N1 = K + N2 = N - K + CALL DLAMRG( N1, N2, D, 1, -1, INDXQ ) + ELSE + DO 60 I = 1, N + INDXQ( I ) = I + 60 CONTINUE + END IF +* +* Persist the parent boundary rows in the current physical D order. +* + IF( WANTQ ) THEN + DO 65 J = 1, N + BLO( J ) = WORK( IQ+2*( J-1 ) ) + BHI( J ) = WORK( IQ+2*( J-1 )+1 ) + 65 CONTINUE + END IF +* + 70 CONTINUE + RETURN +* +* End of DLAED7_BR +* + END +* +*> \brief \b DLAED8_BR selected-row variant of DLAED8. +* +* DLAED8_BR mirrors DLAED8 deflation and bookkeeping but applies the +* column rotations/permutations to NBR selected rows rather than to a +* full eigenvector matrix. It is an internal helper for DLAED7_BR. +* + SUBROUTINE DLAED8_BR( NBR, K, N, D, Q, LDQ, INDXQ, RHO, + $ CUTPNT, Z, DLAMBDA, Q2, LDQ2, W, PERM, + $ GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, + $ INFO ) +* +* .. Scalar Arguments .. + INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, NBR + DOUBLE PRECISION RHO +* .. +* .. Array Arguments .. + INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), + $ INDXQ( * ), PERM( * ) + DOUBLE PRECISION D( * ), DLAMBDA( * ), GIVNUM( 2, * ), + $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT + PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, + $ ONE = 1.0D0, TWO = 2.0D0, EIGHT = 8.0D0 ) +* .. +* .. Local Scalars .. + INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2 + DOUBLE PRECISION C, EPS, S, T, TAU, TOL +* .. +* .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH, DLAPY2 + EXTERNAL IDAMAX, DLAMCH, DLAPY2 +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( NBR.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDQ.LT.MAX( 1, NBR ) ) THEN + INFO = -6 + ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN + INFO = -9 + ELSE IF( LDQ2.LT.MAX( 1, NBR ) ) THEN + INFO = -13 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAED8_BR', -INFO ) + RETURN + END IF +* + GIVPTR = 0 + IF( N.EQ.0 ) + $ RETURN +* + N1 = CUTPNT + N2 = N - N1 + N1P1 = N1 + 1 +* + IF( RHO.LT.ZERO ) THEN + CALL DSCAL( N2, MONE, Z( N1P1 ), 1 ) + END IF +* + T = ONE / SQRT( TWO ) + DO 80 J = 1, N + INDX( J ) = J + 80 CONTINUE + CALL DSCAL( N, T, Z, 1 ) + RHO = ABS( TWO*RHO ) +* +* Sort the eigenvalues into increasing order. +* + DO 90 I = CUTPNT + 1, N + INDXQ( I ) = INDXQ( I ) + CUTPNT + 90 CONTINUE + DO 100 I = 1, N + DLAMBDA( I ) = D( INDXQ( I ) ) + W( I ) = Z( INDXQ( I ) ) + 100 CONTINUE + CALL DLAMRG( N1, N2, DLAMBDA, 1, 1, INDX ) + DO 110 I = 1, N + D( I ) = DLAMBDA( INDX( I ) ) + Z( I ) = W( INDX( I ) ) + 110 CONTINUE +* + IMAX = IDAMAX( N, Z, 1 ) + JMAX = IDAMAX( N, D, 1 ) + EPS = DLAMCH( 'Epsilon' ) + TOL = EIGHT*EPS*ABS( D( JMAX ) ) +* + IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN + K = 0 + IF( NBR.EQ.0 ) THEN + DO 115 J = 1, N + PERM( J ) = INDXQ( INDX( J ) ) + 115 CONTINUE + ELSE + DO 120 J = 1, N + PERM( J ) = INDXQ( INDX( J ) ) + CALL DCOPY( NBR, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) + 120 CONTINUE + CALL DLACPY( 'A', NBR, N, Q2, LDQ2, Q, LDQ ) + END IF + RETURN + END IF +* + K = 0 + K2 = N + 1 + DO 130 J = 1, N + IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN + K2 = K2 - 1 + INDXP( K2 ) = J + IF( J.EQ.N ) + $ GO TO 170 + ELSE + JLAM = J + GO TO 140 + END IF + 130 CONTINUE + 140 CONTINUE + J = J + 1 + IF( J.GT.N ) + $ GO TO 160 + IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN + K2 = K2 - 1 + INDXP( K2 ) = J + ELSE + S = Z( JLAM ) + C = Z( J ) + TAU = DLAPY2( C, S ) + T = D( J ) - D( JLAM ) + C = C / TAU + S = -S / TAU + IF( ABS( T*C*S ).LE.TOL ) THEN + Z( J ) = TAU + Z( JLAM ) = ZERO +* + GIVPTR = GIVPTR + 1 + GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) ) + GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) ) + GIVNUM( 1, GIVPTR ) = C + GIVNUM( 2, GIVPTR ) = S + IF( NBR.NE.0 ) + $ CALL DROT( NBR, Q( 1, INDXQ( INDX( JLAM ) ) ), 1, + $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S ) + T = D( JLAM )*C*C + D( J )*S*S + D( J ) = D( JLAM )*S*S + D( J )*C*C + D( JLAM ) = T + K2 = K2 - 1 + I = 1 + 150 CONTINUE + IF( K2+I.LE.N ) THEN + IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN + INDXP( K2+I-1 ) = INDXP( K2+I ) + INDXP( K2+I ) = JLAM + I = I + 1 + GO TO 150 + ELSE + INDXP( K2+I-1 ) = JLAM + END IF + ELSE + INDXP( K2+I-1 ) = JLAM + END IF + JLAM = J + ELSE + K = K + 1 + W( K ) = Z( JLAM ) + DLAMBDA( K ) = D( JLAM ) + INDXP( K ) = JLAM + JLAM = J + END IF + END IF + GO TO 140 + 160 CONTINUE +* + K = K + 1 + W( K ) = Z( JLAM ) + DLAMBDA( K ) = D( JLAM ) + INDXP( K ) = JLAM +* + 170 CONTINUE + DO 180 J = 1, N + JP = INDXP( J ) + DLAMBDA( J ) = D( JP ) + PERM( J ) = INDXQ( INDX( JP ) ) + IF( NBR.NE.0 ) + $ CALL DCOPY( NBR, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) + 180 CONTINUE +* + IF( K.LT.N ) THEN + CALL DCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) + IF( NBR.NE.0 ) + $ CALL DLACPY( 'A', NBR, N-K, Q2( 1, K+1 ), LDQ2, + $ Q( 1, K+1 ), LDQ ) + END IF +* + RETURN +* +* End of DLAED8_BR +* + END diff --git a/lapack-netlib/SRC/dlaed9_br.f b/lapack-netlib/SRC/dlaed9_br.f new file mode 100644 index 0000000000..38d153f6c2 --- /dev/null +++ b/lapack-netlib/SRC/dlaed9_br.f @@ -0,0 +1,253 @@ +*> \brief \b DLAED9_BR used by DSTEBR. +* +* DLAED9_BR finds the roots of the secular equation for one +* boundary-row DC merge and directly updates a small set of selected +* rows. It mirrors DLAED9's root solve and eigenvector normalization, +* but avoids materializing the full K-by-K secular eigenvector block. +* DLAED4_BR returns high-relative-accuracy DELTA values for each root; +* consume those DELTA values immediately to form WPROD, following the +* same numerical idea as LAPACK's SVD divide-and-conquer helpers. +* +* Serial Reference-LAPACK form: no OpenMP or thread-private accumulators. +* + SUBROUTINE DLAED9_BR( K, NBR, D, QIN, LDQIN, QOUT, LDQOUT, + $ RHO, DLAMBDA, W, DELTA, WPROD, WSGN, + $ TAU, NEAR, INFO ) +* +* .. Scalar Arguments .. + INTEGER INFO, K, LDQIN, LDQOUT, NBR + DOUBLE PRECISION RHO +* .. +* .. Array Arguments .. + DOUBLE PRECISION D( * ), DELTA( * ), DLAMBDA( * ), + $ NEAR( 2, * ), QIN( LDQIN, * ), + $ QOUT( LDQOUT, * ), TAU( * ), W( * ), + $ WPROD( * ), WSGN( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +* .. +* .. Local Scalars .. + INTEGER I, IORG, J, N1, N2, R + DOUBLE PRECISION ACC1, ACC2, DLT, ROWSUM, TEMP +* .. +* .. External Subroutines .. + EXTERNAL DLAED4, DLAED4_BR, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SIGN, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( K.LT.0 ) THEN + INFO = -1 + ELSE IF( NBR.LT.0 ) THEN + INFO = -2 + ELSE IF( LDQIN.LT.MAX( 1, NBR ) ) THEN + INFO = -5 + ELSE IF( LDQOUT.LT.MAX( 1, NBR ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAED9_BR', -INFO ) + RETURN + END IF +* + IF( K.EQ.0 ) + $ RETURN +* +* Root-only final merge: no parent boundary rows will be observed by +* a higher merge level, so updated W and secular vectors are not +* needed. +* + IF( NBR.EQ.0 ) THEN + DO 5 J = 1, K + CALL DLAED4( K, J, DLAMBDA, W, DELTA, RHO, D( J ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 170 + 5 CONTINUE + GO TO 170 + END IF +* +* For K = 1 or 2, DLAED4/DLAED5 returns normalized eigenvectors +* directly in DELTA. Keep this path identical to DLAED9's +* special-case copy to S. +* + IF( K.EQ.1 .OR. K.EQ.2 ) THEN + DO 30 J = 1, K + CALL DLAED4( K, J, DLAMBDA, W, DELTA, RHO, D( J ), + $ INFO ) + IF( INFO.NE.0 ) + $ GO TO 170 + IF( NBR.EQ.2 ) THEN + ACC1 = ZERO + ACC2 = ZERO + DO 10 I = 1, K + ACC1 = ACC1 + QIN( 1, I )*DELTA( I ) + ACC2 = ACC2 + QIN( 2, I )*DELTA( I ) + 10 CONTINUE + QOUT( 1, J ) = ACC1 + QOUT( 2, J ) = ACC2 + ELSE + DO 20 R = 1, NBR + ROWSUM = ZERO + DO 15 I = 1, K + ROWSUM = ROWSUM + QIN( R, I )*DELTA( I ) + 15 CONTINUE + QOUT( R, J ) = ROWSUM + 20 CONTINUE + END IF + 30 CONTINUE + GO TO 170 + END IF +* +* Preserve the update vector, solve all secular roots, then accumulate +* the product needed to reconstruct DLAED9's updated W. The +* high-relative-accuracy DELTA values produced by DLAED4_BR are +* consumed immediately rather than reconstructed later from the +* compact origin/tau representation. +* + DO 40 I = 1, K + WSGN( I ) = W( I ) + 40 CONTINUE +* + DO 41 I = 1, K + WPROD( I ) = ONE + 41 CONTINUE +* + DO 44 J = 1, K + CALL DLAED4_BR( K, J, DLAMBDA, WSGN, DELTA, RHO, D( J ), + $ IORG, TAU( J ), NEAR( 1, J ), + $ NEAR( 2, J ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 170 +* + WPROD( J ) = WPROD( J )*DELTA( J ) + DO 42 I = 1, J - 1 + WPROD( I ) = WPROD( I )* + $ ( DELTA( I ) / + $ ( DLAMBDA( I )-DLAMBDA( J ) ) ) + 42 CONTINUE + DO 43 I = J + 1, K + WPROD( I ) = WPROD( I )* + $ ( DELTA( I ) / + $ ( DLAMBDA( I )-DLAMBDA( J ) ) ) + 43 CONTINUE + 44 CONTINUE +* + DO 80 I = 1, K + W( I ) = SIGN( SQRT( -WPROD( I ) ), WSGN( I ) ) + 80 CONTINUE +* +* Stream each secular eigenvector column through the selected rows. +* Reconstruct the DELTA column from the compact origin/tau state, +* using cached pole-adjacent values where cancellation is largest. +* + IF( NBR.EQ.2 ) THEN + DO 115 J = 1, K + IF( J.EQ.K ) THEN + IORG = K + N1 = K - 1 + N2 = K + ELSE IF( TAU( J ).GE.ZERO ) THEN + IORG = J + N1 = J + N2 = J + 1 + ELSE + IORG = J + 1 + N1 = J + N2 = J + 1 + END IF + TEMP = ZERO + ACC1 = ZERO + ACC2 = ZERO +* + DO 85 I = 1, N1 - 1 + DLT = W( I ) / + $ ( ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, I )*DLT + ACC2 = ACC2 + QIN( 2, I )*DLT + 85 CONTINUE +* + DLT = W( N1 ) / NEAR( 1, J ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, N1 )*DLT + ACC2 = ACC2 + QIN( 2, N1 )*DLT +* + DO 86 I = N1 + 1, N2 - 1 + DLT = W( I ) / + $ ( ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, I )*DLT + ACC2 = ACC2 + QIN( 2, I )*DLT + 86 CONTINUE +* + DLT = W( N2 ) / NEAR( 2, J ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, N2 )*DLT + ACC2 = ACC2 + QIN( 2, N2 )*DLT +* + DO 87 I = N2 + 1, K + DLT = W( I ) / + $ ( ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, I )*DLT + ACC2 = ACC2 + QIN( 2, I )*DLT + 87 CONTINUE +* + TEMP = SQRT( TEMP ) + QOUT( 1, J ) = ACC1 / TEMP + QOUT( 2, J ) = ACC2 / TEMP + 115 CONTINUE + GO TO 170 + END IF +* + DO 130 J = 1, K + IF( J.EQ.K ) THEN + IORG = K + ELSE IF( TAU( J ).GE.ZERO ) THEN + IORG = J + ELSE + IORG = J + 1 + END IF + TEMP = ZERO + DO 90 R = 1, NBR + QOUT( R, J ) = ZERO + 90 CONTINUE + DO 100 I = 1, K + DLT = ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) + IF( J.EQ.K ) THEN + IF( I.EQ.K-1 ) + $ DLT = NEAR( 1, J ) + IF( I.EQ.K ) + $ DLT = NEAR( 2, J ) + ELSE + IF( I.EQ.J ) + $ DLT = NEAR( 1, J ) + IF( I.EQ.J+1 ) + $ DLT = NEAR( 2, J ) + END IF + DLT = W( I ) / DLT + TEMP = TEMP + DLT*DLT + DO 95 R = 1, NBR + QOUT( R, J ) = QOUT( R, J ) + QIN( R, I )*DLT + 95 CONTINUE + 100 CONTINUE + TEMP = SQRT( TEMP ) + DO 120 R = 1, NBR + QOUT( R, J ) = QOUT( R, J ) / TEMP + 120 CONTINUE + 130 CONTINUE +* + 170 CONTINUE + RETURN +* +* End of DLAED9_BR +* + END diff --git a/lapack-netlib/SRC/dstebr.f b/lapack-netlib/SRC/dstebr.f new file mode 100644 index 0000000000..dfa7f8ef5e --- /dev/null +++ b/lapack-netlib/SRC/dstebr.f @@ -0,0 +1,274 @@ +*> \brief \b DSTEBR +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DSTEBR( N, D, E, WORK, LWORK, IWORK, LIWORK, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LIWORK, LWORK, N +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION D( * ), E( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSTEBR computes all eigenvalues of a symmetric tridiagonal matrix +*> using a values-only divide-and-conquer method that propagates only +*> the first and last boundary rows of each local eigenvector block +*> (STEBR = STymmetric tridiagonal Eigenvalues, Boundary-Row DC). +*> +*> Unlike DSTEDC with COMPZ = 'N', DSTEBR does not store full secular +*> eigenvector blocks across merge levels for later replay. It is a +*> distinct routine and does not replace DSTEDC or DSTERF. +*> +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix. N >= 0. +*> \endverbatim +*> +*> \param[in,out] D +*> \verbatim +*> D is DOUBLE PRECISION array, dimension (N) +*> On entry, the n diagonal elements of the tridiagonal matrix. +*> On exit, if INFO = 0, the eigenvalues in ascending order. +*> \endverbatim +*> +*> \param[in,out] E +*> \verbatim +*> E is DOUBLE PRECISION array, dimension (N-1) +*> On entry, the (n-1) subdiagonal elements of the tridiagonal +*> matrix. +*> On exit, E has been destroyed. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK. +*> If N <= SMLSIZ, where SMLSIZ is returned by ILAENV and is +*> typically about 25, LWORK must be at least 1. +*> If N > SMLSIZ, LWORK must be at least +*> MAX( 16*N, 2*N + SMLSIZ**2 + 2*SMLSIZ - 2 ). +*> +*> If LWORK = -1, then a workspace query is assumed; the +*> routine only calculates the optimal sizes of the WORK and +*> IWORK arrays, returns these values as the first entries of +*> the WORK and IWORK arrays, and no error message related to +*> LWORK or LIWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) +*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. +*> \endverbatim +*> +*> \param[in] LIWORK +*> \verbatim +*> LIWORK is INTEGER +*> The dimension of the array IWORK. +*> If N <= SMLSIZ, LIWORK must be at least 1. +*> If N > SMLSIZ, LIWORK must be at least 7*N. +*> +*> If LIWORK = -1, then a workspace query is assumed; the +*> routine only calculates the optimal sizes of the WORK and +*> IWORK arrays, returns these values as the first entries of +*> the WORK and IWORK arrays, and no error message related to +*> LWORK or LIWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: The algorithm failed to compute an eigenvalue while +*> working on a submatrix lying in rows and columns +*> INFO/(N+1) through MOD(INFO,N+1). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup auxOTHERcomputational +* +* ===================================================================== + SUBROUTINE DSTEBR( N, D, E, WORK, LWORK, IWORK, LIWORK, INFO ) +* +* -- LAPACK computational routine -- +* +* .. Scalar Arguments .. + INTEGER INFO, LIWORK, LWORK, N +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION D( * ), E( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +* .. +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER FINISH, LIWMIN, LWMIN, M, SMLSIZ, START + DOUBLE PRECISION EPS, ORGNRM, TINY +* .. +* .. External Functions .. + INTEGER ILAENV + DOUBLE PRECISION DLAMCH, DLANST + EXTERNAL ILAENV, DLAMCH, DLANST +* .. +* .. External Subroutines .. + EXTERNAL DLAED0_BR, DLASCL, DLASRT, DSTERF, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MOD, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) +* + IF( N.LT.0 ) THEN + INFO = -1 + END IF +* + IF( INFO.EQ.0 ) THEN + SMLSIZ = ILAENV( 9, 'DSTEBR', ' ', 0, 0, 0, 0 ) + IF( N.LE.SMLSIZ ) THEN + LWMIN = 1 + LIWMIN = 1 + ELSE +* +* BLO/BHI: 2*N; merge scratch: 14*N (see DLAED7_BR). +* Leaf solves need 2*N + SMLSIZ**2 + 2*SMLSIZ - 2: the +* SMLSIZ-by-SMLSIZ leaf eigenvector block plus the +* DSTEQR work array. +* + LWMIN = MAX( 16*N, 2*N + SMLSIZ*SMLSIZ + 2*SMLSIZ - 2 ) + LIWMIN = MAX( 1, 7*N ) + END IF + WORK( 1 ) = LWMIN + IWORK( 1 ) = LIWMIN +* + IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN + INFO = -5 + ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSTEBR', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN + IF( N.EQ.1 ) + $ RETURN +* + IF( N.LE.SMLSIZ ) THEN + CALL DSTERF( N, D, E, INFO ) + GO TO 50 + END IF +* + ORGNRM = DLANST( 'M', N, D, E ) + IF( ORGNRM.LE.ZERO ) + $ GO TO 50 +* + EPS = DLAMCH( 'Epsilon' ) + START = 1 +* + 10 CONTINUE + IF( START.LE.N ) THEN + FINISH = START + 20 CONTINUE + IF( FINISH.LT.N ) THEN + TINY = EPS*SQRT( ABS( D( FINISH ) ) )* + $ SQRT( ABS( D( FINISH+1 ) ) ) + IF( ABS( E( FINISH ) ).GT.TINY ) THEN + FINISH = FINISH + 1 + GO TO 20 + END IF + END IF +* + M = FINISH - START + 1 + IF( M.EQ.1 ) THEN + START = FINISH + 1 + GO TO 10 + END IF +* + IF( M.GT.SMLSIZ ) THEN + ORGNRM = DLANST( 'M', M, D( START ), E( START ) ) + CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), + $ M, INFO ) + CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, + $ E( START ), M-1, INFO ) +* + CALL DLAED0_BR( M, D( START ), E( START ), WORK, IWORK, + $ INFO ) + IF( INFO.NE.0 ) THEN + INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + + $ MOD( INFO, ( M+1 ) ) + START - 1 + GO TO 50 + END IF +* + CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), + $ M, INFO ) + ELSE + CALL DSTERF( M, D( START ), E( START ), INFO ) + IF( INFO.NE.0 ) THEN + INFO = START*( N+1 ) + FINISH + GO TO 50 + END IF + END IF +* + START = FINISH + 1 + GO TO 10 + END IF +* + CALL DLASRT( 'I', N, D, INFO ) +* + 50 CONTINUE + WORK( 1 ) = LWMIN + IWORK( 1 ) = LIWMIN + RETURN +* +* End of DSTEBR +* + END diff --git a/lapack-netlib/SRC/slaed0_br.f b/lapack-netlib/SRC/slaed0_br.f new file mode 100644 index 0000000000..a25f69a6a1 --- /dev/null +++ b/lapack-netlib/SRC/slaed0_br.f @@ -0,0 +1,233 @@ +*> \brief \b SLAED0_BR used by SSTEBR. +* +* SLAED0_BR computes eigenvalues of a symmetric tridiagonal matrix with +* the boundary-row divide-and-conquer method used by SSTEBR. Unlike +* SLAED0(ICOMPQ=0), it keeps only first/last boundary rows for each +* active block across merge levels instead of saving complete +* eigenvector blocks. +* +* Serial Reference-LAPACK form: no OpenMP, environment-variable +* scheduling, or thread-scaled workspace. +* + SUBROUTINE SLAED0_BR( N, D, E, WORK, IWORK, INFO ) +* +* .. Scalar Arguments .. + INTEGER INFO, N +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL D( * ), E( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + INTEGER WRKSTR + PARAMETER ( WRKSTR = 14 ) +* .. +* .. Local Scalars .. + LOGICAL WANTQ + INTEGER I, IBHI, IBLO, IIWRK, INDXQ, INFOM, ISCR, + $ IWBASE, J, K, KDEFL, MATSIZ, MERGE, MSD2, + $ NMERGE, QBASE, SMLSIZ, SMM1, SPM1, SUBMAT, + $ SUBPBS, WBASE, WRKBASE +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SLAED7_BR, SSTEQR_BOUNDARY_BR, XERBLA +* .. +* .. External Functions .. + INTEGER ILAENV + EXTERNAL ILAENV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( N.LT.0 ) THEN + INFO = -1 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLAED0_BR', -INFO ) + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN +* + SMLSIZ = ILAENV( 9, 'SLAED0', ' ', 0, 0, 0, 0 ) +* +* Determine the size and placement of leaf submatrices. +* + IWORK( 1 ) = N + SUBPBS = 1 + 10 CONTINUE + IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN + DO 20 J = SUBPBS, 1, -1 + IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 + IWORK( 2*J-1 ) = IWORK( J ) / 2 + 20 CONTINUE + SUBPBS = 2*SUBPBS + GO TO 10 + END IF + DO 30 J = 2, SUBPBS + IWORK( J ) = IWORK( J ) + IWORK( J-1 ) + 30 CONTINUE +* +* Apply rank-1 cuts. +* + SPM1 = SUBPBS - 1 + DO 40 I = 1, SPM1 + SUBMAT = IWORK( I ) + 1 + SMM1 = SUBMAT - 1 + D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) + D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) + 40 CONTINUE +* +* Workspace layout: +* WORK(IBLO:IBLO+N-1) first boundary rows +* WORK(IBHI:IBHI+N-1) last boundary rows +* WORK(ISCR:) leaf scratch and per-merge scratch +* +* Integer workspace layout: +* IWORK(1:N) split endpoints +* IWORK(INDXQ:INDXQ+N-1) sorting permutations +* IWORK(IIWRK:) SLAED7_BR scratch +* + IBLO = 1 + IBHI = IBLO + N + ISCR = IBHI + N + INDXQ = N + 1 + IIWRK = INDXQ + N +* +* Solve leaf eigenproblems and initialize boundary row invariants. +* + DO 70 I = 1, SUBPBS + QBASE = ISCR + WBASE = QBASE + SMLSIZ*SMLSIZ + IF( I.EQ.1 ) THEN + SUBMAT = 1 + MATSIZ = IWORK( 1 ) + ELSE + SUBMAT = IWORK( I-1 ) + 1 + MATSIZ = IWORK( I ) - IWORK( I-1 ) + END IF +* + CALL SSTEQR_BOUNDARY_BR( MATSIZ, D( SUBMAT ), E( SUBMAT ), + $ WORK( IBLO+SUBMAT-1 ), + $ WORK( IBHI+SUBMAT-1 ), + $ WORK( QBASE ), SMLSIZ, WORK( WBASE ), + $ INFO ) + IF( INFO.NE.0 ) + $ GO TO 130 + K = 1 + DO 65 J = SUBMAT, IWORK( I ) + IWORK( INDXQ+J-1 ) = K + K = K + 1 + 65 CONTINUE + 70 CONTINUE +* +* Merge adjacent eigensystems. Boundary rows are updated in place +* for the merged physical D order. +* + 80 CONTINUE + IF( SUBPBS.GT.1 ) THEN + WANTQ = SUBPBS.GT.2 + NMERGE = SUBPBS / 2 + DO 90 MERGE = 1, NMERGE + CALL SLAED0_MERGE_PAIR( MERGE, NMERGE, IWORK, WANTQ, D, E, + $ WORK, IBLO, IBHI, ISCR, IIWRK, INDXQ, SUBMAT, MATSIZ, + $ MSD2, WRKBASE, IWBASE, KDEFL, INFOM, WRKSTR ) + IF( INFOM.NE.0 ) THEN + INFO = INFOM + GO TO 130 + END IF + 90 CONTINUE +* + DO 95 MERGE = 1, NMERGE + IWORK( MERGE ) = IWORK( 2*MERGE ) + 95 CONTINUE + SUBPBS = SUBPBS / 2 + GO TO 80 + END IF +* +* Reorder eigenvalues from the final physical merge order. +* + DO 100 I = 1, N + J = IWORK( INDXQ+I-1 ) + WORK( ISCR+I-1 ) = D( J ) + 100 CONTINUE + CALL SCOPY( N, WORK( ISCR ), 1, D, 1 ) + GO TO 140 +* + 130 CONTINUE + INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 +* + 140 CONTINUE + RETURN +* +* End of SLAED0_BR +* + END +* +*> \brief \b SLAED0_MERGE_PAIR executes one merge in the DC tree. +* + SUBROUTINE SLAED0_MERGE_PAIR( MERGE, NMERGE, IWORK, WANTQ, D, E, + $ WORK, IBLO, IBHI, ISCR, IIWRK, + $ INDXQ, SUBMAT, MATSIZ, MSD2, + $ WRKBASE, IWBASE, KDEFL, INFO, + $ WRKSTR ) +* + INTEGER IBHI, IBLO, IIWRK, INDXQ, INFO, ISCR, IWBASE, + $ KDEFL, MATSIZ, MERGE, MSD2, NMERGE, RIGHT, + $ SUBMAT, WRKBASE, WRKSTR + LOGICAL WANTQ + INTEGER IWORK( * ) + REAL D( * ), E( * ), WORK( * ) + EXTERNAL SLAED7_BR +* + INTEGER LEFT +* + LEFT = 2*MERGE - 1 + RIGHT = LEFT + 1 + IF( MERGE.EQ.1 ) THEN + SUBMAT = 1 + MATSIZ = IWORK( RIGHT ) + MSD2 = IWORK( 1 ) + ELSE + SUBMAT = IWORK( LEFT-1 ) + 1 + MATSIZ = IWORK( RIGHT ) - IWORK( LEFT-1 ) + MSD2 = IWORK( LEFT ) - IWORK( LEFT-1 ) + END IF +* + WRKBASE = ISCR + WRKSTR*( SUBMAT-1 ) + IWBASE = IIWRK + 5*( SUBMAT-1 ) +* + CALL SLAED7_BR( WANTQ, MATSIZ, MSD2, D( SUBMAT ), + $ IWORK( INDXQ+SUBMAT-1 ), E( SUBMAT+MSD2-1 ), + $ WORK( IBLO+SUBMAT-1 ), WORK( IBHI+SUBMAT-1 ), + $ WORK( WRKBASE ), IWORK( IWBASE ), KDEFL, INFO ) + RETURN + END +* +*> \brief \b SSTEQR_BOUNDARY_BR solves one leaf and keeps boundary rows. +* + SUBROUTINE SSTEQR_BOUNDARY_BR( N, D, E, BLO, BHI, Q, LDQ, WORK, + $ INFO ) +* + INTEGER INFO, J, LDQ, N + REAL BLO( * ), BHI( * ), D( * ), E( * ), Q( * ), + $ WORK( * ) + EXTERNAL SSTEQR +* + CALL SSTEQR( 'I', N, D, E, Q, LDQ, WORK, INFO ) + IF( INFO.NE.0 ) + $ RETURN +* + DO 10 J = 1, N + BLO( J ) = Q( 1 + ( J - 1 )*LDQ ) + BHI( J ) = Q( N + ( J - 1 )*LDQ ) + 10 CONTINUE + RETURN + END diff --git a/lapack-netlib/SRC/slaed4_br.f b/lapack-netlib/SRC/slaed4_br.f new file mode 100644 index 0000000000..a74ff427ac --- /dev/null +++ b/lapack-netlib/SRC/slaed4_br.f @@ -0,0 +1,109 @@ +*> \brief \b SLAED4_BR used by SSTEBR. +* +* SLAED4_BR computes one secular root with SLAED4 and exposes a compact +* representation of the returned high-relative-accuracy DELTA vector. +* For N > 2, SLAED4's DELTA entries are all relative to either D(I) or +* D(I+1), except for the final root where the origin is D(N). ORG records +* that origin index and TAU satisfies +* +* lambda = D(ORG) + TAU, +* DELTA(j) = D(j) - D(ORG) - TAU. +* +* NEAR1/NEAR2 cache the two pole-adjacent DELTA entries that are most +* sensitive to cancellation when reconstructed from D and TAU. +* + SUBROUTINE SLAED4_BR( N, I, D, Z, DELTA, RHO, DLAM, ORG, + $ TAU, NEAR1, NEAR2, INFO ) +* +* .. Scalar Arguments .. + INTEGER I, INFO, N, ORG + REAL DLAM, NEAR1, NEAR2, RHO, TAU +* .. +* .. Array Arguments .. + REAL D( * ), DELTA( * ), Z( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL HALF, ZERO + PARAMETER ( HALF = 0.5E0, ZERO = 0.0E0 ) +* .. +* .. Local Scalars .. + REAL MID +* .. +* .. External Subroutines .. + EXTERNAL SLAED4, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS +* .. +* .. Executable Statements .. +* + INFO = 0 + ORG = 1 + TAU = ZERO + NEAR1 = ZERO + NEAR2 = ZERO +* + IF( N.LT.1 ) THEN + INFO = -1 + ELSE IF( I.LT.1 .OR. I.GT.N ) THEN + INFO = -2 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLAED4_BR', -INFO ) + RETURN + END IF +* + CALL SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) + IF( INFO.NE.0 ) + $ RETURN +* + IF( N.EQ.1 ) THEN + ORG = 1 + TAU = DLAM - D( 1 ) + NEAR1 = DELTA( 1 ) + NEAR2 = ZERO + RETURN + END IF +* + IF( N.EQ.2 ) THEN + ORG = 1 + TAU = DLAM - D( 1 ) + NEAR1 = DELTA( 1 ) + NEAR2 = DELTA( 2 ) + RETURN + END IF +* + IF( I.EQ.N ) THEN + ORG = N + ELSE + MID = HALF*( D( I )+D( I+1 ) ) + IF( DLAM.LT.MID ) THEN + ORG = I + ELSE + ORG = I + 1 + END IF + END IF +* +* Avoid reconstructing TAU from DLAM-D(ORG). SLAED4 already +* returned DELTA(ORG) in its shifted coordinate, so using +* TAU=-DELTA(ORG) preserves that final rounded denominator and +* avoids an extra add/subtract rounding through DLAM. +* + TAU = -DELTA( ORG ) +* + IF( I.EQ.N ) THEN + NEAR1 = DELTA( N-1 ) + NEAR2 = DELTA( N ) + ELSE + NEAR1 = DELTA( I ) + NEAR2 = DELTA( I+1 ) + END IF +* + RETURN +* +* End of SLAED4_BR +* + END diff --git a/lapack-netlib/SRC/slaed7_br.f b/lapack-netlib/SRC/slaed7_br.f new file mode 100644 index 0000000000..2c6cb28871 --- /dev/null +++ b/lapack-netlib/SRC/slaed7_br.f @@ -0,0 +1,376 @@ +*> \brief \b SLAED7_BR used by SSTEBR. +* +* SLAED7_BR performs one divide-and-conquer merge while propagating +* only the first and last rows of the current local eigenvector block. +* It is called by SLAED0_BR and is not itself a public LAPACK API. +* +* Serial Reference-LAPACK form: no OpenMP or thread-scaled workspace. +* + SUBROUTINE SLAED7_BR( WANTQ, N, CUTPNT, D, INDXQ, RHO, BLO, BHI, + $ WORK, IWORK, KDEFL, INFO ) +* +* .. Scalar Arguments .. + LOGICAL WANTQ + INTEGER CUTPNT, INFO, KDEFL, N + REAL RHO +* .. +* .. Array Arguments .. + INTEGER INDXQ( * ), IWORK( * ) + REAL BLO( * ), BHI( * ), D( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0E0 ) +* .. +* .. Local Scalars .. + INTEGER GIVPTR, IDLMDA, IGIVCL, IGIVNM, I, IINDX, + $ IINDXP, INEAR, IPERM, IQ, IQ2, ISCR, ITAU, IW, + $ IWSGN, IZ, J, K, N1, N2, NBR +* .. +* .. External Subroutines .. + EXTERNAL SLAED8_BR, SLAED9_BR, SLAMRG, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* + INFO = 0 + KDEFL = 0 + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN + INFO = -2 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLAED7_BR', -INFO ) + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN +* +* Workspace layout (14*N doubles): +* +* Z, DLAMBDA, W : N each +* QBR, Q2BR, GIVNUM/reused WPROD : 2*N each +* DELTA, WSGN, TAU, NEAR : N, N, N, 2*N +* + IZ = 1 + IDLMDA = IZ + N + IW = IDLMDA + N + IQ = IW + N + IQ2 = IQ + 2*N + IGIVNM = IQ2 + 2*N + ISCR = IGIVNM + 2*N + IWSGN = ISCR + N + ITAU = IWSGN + N + INEAR = ITAU + N +* +* Integer workspace layout. +* + IINDXP = 1 + IINDX = IINDXP + N + IPERM = IINDX + N + IGIVCL = IPERM + N +* + IF( WANTQ ) THEN + NBR = 2 + ELSE + NBR = 0 + END IF +* +* Z is the rank-one update vector. QBR has two selected rows of +* blockdiag(Q_left,Q_right): parent first boundary and parent last +* boundary, represented in the child eigenbases. The root merge has +* no parent, so it only needs Z for the secular eigenvalues. +* + IF( WANTQ ) THEN + DO 20 J = 1, CUTPNT + WORK( IZ+J-1 ) = BHI( J ) + WORK( IQ+2*( J-1 ) ) = BLO( J ) + WORK( IQ+2*( J-1 )+1 ) = ZERO + 20 CONTINUE + DO 30 J = CUTPNT + 1, N + WORK( IZ+J-1 ) = BLO( J ) + WORK( IQ+2*( J-1 ) ) = ZERO + WORK( IQ+2*( J-1 )+1 ) = BHI( J ) + 30 CONTINUE + ELSE + DO 40 J = 1, CUTPNT + WORK( IZ+J-1 ) = BHI( J ) + 40 CONTINUE + DO 50 J = CUTPNT + 1, N + WORK( IZ+J-1 ) = BLO( J ) + 50 CONTINUE + END IF +* +* Sort and deflate, applying the same Givens/permutation actions to +* the selected boundary rows when a parent merge will need them. +* +* + CALL SLAED8_BR( NBR, K, N, D, WORK( IQ ), 2, INDXQ, RHO, + $ CUTPNT, WORK( IZ ), WORK( IDLMDA ), + $ WORK( IQ2 ), 2, WORK( IW ), IWORK( IPERM ), + $ GIVPTR, IWORK( IGIVCL ), WORK( IGIVNM ), + $ IWORK( IINDXP ), IWORK( IINDX ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 70 + KDEFL = K +* +* Solve the secular equation and directly update the selected rows. +* WORK(ISCR) holds one SLAED4 delta vector at a time; no K-by-K +* delta or eigenvector block is materialized. +* + IF( K.NE.0 ) THEN + CALL SLAED9_BR( K, NBR, D, WORK( IQ2 ), 2, WORK( IQ ), 2, + $ RHO, WORK( IDLMDA ), WORK( IW ), + $ WORK( ISCR ), WORK( IGIVNM ), + $ WORK( IWSGN ), WORK( ITAU ), + $ WORK( INEAR ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 70 +* +* Prepare INDXQ sorting permutation. +* + N1 = K + N2 = N - K + CALL SLAMRG( N1, N2, D, 1, -1, INDXQ ) + ELSE + DO 60 I = 1, N + INDXQ( I ) = I + 60 CONTINUE + END IF +* +* Persist the parent boundary rows in the current physical D order. +* + IF( WANTQ ) THEN + DO 65 J = 1, N + BLO( J ) = WORK( IQ+2*( J-1 ) ) + BHI( J ) = WORK( IQ+2*( J-1 )+1 ) + 65 CONTINUE + END IF +* + 70 CONTINUE + RETURN +* +* End of SLAED7_BR +* + END +* +*> \brief \b SLAED8_BR selected-row variant of SLAED8. +* +* SLAED8_BR mirrors SLAED8 deflation and bookkeeping but applies the +* column rotations/permutations to NBR selected rows rather than to a +* full eigenvector matrix. It is an internal helper for SLAED7_BR. +* + SUBROUTINE SLAED8_BR( NBR, K, N, D, Q, LDQ, INDXQ, RHO, + $ CUTPNT, Z, DLAMBDA, Q2, LDQ2, W, PERM, + $ GIVPTR, GIVCOL, GIVNUM, INDXP, INDX, + $ INFO ) +* +* .. Scalar Arguments .. + INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, NBR + REAL RHO +* .. +* .. Array Arguments .. + INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), + $ INDXQ( * ), PERM( * ) + REAL D( * ), DLAMBDA( * ), GIVNUM( 2, * ), + $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL MONE, ZERO, ONE, TWO, EIGHT + PARAMETER ( MONE = -1.0E0, ZERO = 0.0E0, + $ ONE = 1.0E0, TWO = 2.0E0, EIGHT = 8.0E0 ) +* .. +* .. Local Scalars .. + INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2 + REAL C, EPS, S, T, TAU, TOL +* .. +* .. External Functions .. + INTEGER ISAMAX + REAL SLAMCH, SLAPY2 + EXTERNAL ISAMAX, SLAMCH, SLAPY2 +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SLACPY, SLAMRG, SROT, SSCAL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( NBR.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDQ.LT.MAX( 1, NBR ) ) THEN + INFO = -6 + ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN + INFO = -9 + ELSE IF( LDQ2.LT.MAX( 1, NBR ) ) THEN + INFO = -13 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLAED8_BR', -INFO ) + RETURN + END IF +* + GIVPTR = 0 + IF( N.EQ.0 ) + $ RETURN +* + N1 = CUTPNT + N2 = N - N1 + N1P1 = N1 + 1 +* + IF( RHO.LT.ZERO ) THEN + CALL SSCAL( N2, MONE, Z( N1P1 ), 1 ) + END IF +* + T = ONE / SQRT( TWO ) + DO 80 J = 1, N + INDX( J ) = J + 80 CONTINUE + CALL SSCAL( N, T, Z, 1 ) + RHO = ABS( TWO*RHO ) +* +* Sort the eigenvalues into increasing order. +* + DO 90 I = CUTPNT + 1, N + INDXQ( I ) = INDXQ( I ) + CUTPNT + 90 CONTINUE + DO 100 I = 1, N + DLAMBDA( I ) = D( INDXQ( I ) ) + W( I ) = Z( INDXQ( I ) ) + 100 CONTINUE + CALL SLAMRG( N1, N2, DLAMBDA, 1, 1, INDX ) + DO 110 I = 1, N + D( I ) = DLAMBDA( INDX( I ) ) + Z( I ) = W( INDX( I ) ) + 110 CONTINUE +* + IMAX = ISAMAX( N, Z, 1 ) + JMAX = ISAMAX( N, D, 1 ) + EPS = SLAMCH( 'Epsilon' ) + TOL = EIGHT*EPS*ABS( D( JMAX ) ) +* + IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN + K = 0 + IF( NBR.EQ.0 ) THEN + DO 115 J = 1, N + PERM( J ) = INDXQ( INDX( J ) ) + 115 CONTINUE + ELSE + DO 120 J = 1, N + PERM( J ) = INDXQ( INDX( J ) ) + CALL SCOPY( NBR, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) + 120 CONTINUE + CALL SLACPY( 'A', NBR, N, Q2, LDQ2, Q, LDQ ) + END IF + RETURN + END IF +* + K = 0 + K2 = N + 1 + DO 130 J = 1, N + IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN + K2 = K2 - 1 + INDXP( K2 ) = J + IF( J.EQ.N ) + $ GO TO 170 + ELSE + JLAM = J + GO TO 140 + END IF + 130 CONTINUE + 140 CONTINUE + J = J + 1 + IF( J.GT.N ) + $ GO TO 160 + IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN + K2 = K2 - 1 + INDXP( K2 ) = J + ELSE + S = Z( JLAM ) + C = Z( J ) + TAU = SLAPY2( C, S ) + T = D( J ) - D( JLAM ) + C = C / TAU + S = -S / TAU + IF( ABS( T*C*S ).LE.TOL ) THEN + Z( J ) = TAU + Z( JLAM ) = ZERO +* + GIVPTR = GIVPTR + 1 + GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) ) + GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) ) + GIVNUM( 1, GIVPTR ) = C + GIVNUM( 2, GIVPTR ) = S + IF( NBR.NE.0 ) + $ CALL SROT( NBR, Q( 1, INDXQ( INDX( JLAM ) ) ), 1, + $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S ) + T = D( JLAM )*C*C + D( J )*S*S + D( J ) = D( JLAM )*S*S + D( J )*C*C + D( JLAM ) = T + K2 = K2 - 1 + I = 1 + 150 CONTINUE + IF( K2+I.LE.N ) THEN + IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN + INDXP( K2+I-1 ) = INDXP( K2+I ) + INDXP( K2+I ) = JLAM + I = I + 1 + GO TO 150 + ELSE + INDXP( K2+I-1 ) = JLAM + END IF + ELSE + INDXP( K2+I-1 ) = JLAM + END IF + JLAM = J + ELSE + K = K + 1 + W( K ) = Z( JLAM ) + DLAMBDA( K ) = D( JLAM ) + INDXP( K ) = JLAM + JLAM = J + END IF + END IF + GO TO 140 + 160 CONTINUE +* + K = K + 1 + W( K ) = Z( JLAM ) + DLAMBDA( K ) = D( JLAM ) + INDXP( K ) = JLAM +* + 170 CONTINUE + DO 180 J = 1, N + JP = INDXP( J ) + DLAMBDA( J ) = D( JP ) + PERM( J ) = INDXQ( INDX( JP ) ) + IF( NBR.NE.0 ) + $ CALL SCOPY( NBR, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) + 180 CONTINUE +* + IF( K.LT.N ) THEN + CALL SCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) + IF( NBR.NE.0 ) + $ CALL SLACPY( 'A', NBR, N-K, Q2( 1, K+1 ), LDQ2, + $ Q( 1, K+1 ), LDQ ) + END IF +* + RETURN +* +* End of SLAED8_BR +* + END diff --git a/lapack-netlib/SRC/slaed9_br.f b/lapack-netlib/SRC/slaed9_br.f new file mode 100644 index 0000000000..fd62a2979a --- /dev/null +++ b/lapack-netlib/SRC/slaed9_br.f @@ -0,0 +1,253 @@ +*> \brief \b SLAED9_BR used by SSTEBR. +* +* SLAED9_BR finds the roots of the secular equation for one +* boundary-row DC merge and directly updates a small set of selected +* rows. It mirrors SLAED9's root solve and eigenvector normalization, +* but avoids materializing the full K-by-K secular eigenvector block. +* SLAED4_BR returns high-relative-accuracy DELTA values for each root; +* consume those DELTA values immediately to form WPROD, following the +* same numerical idea as LAPACK's SVD divide-and-conquer helpers. +* +* Serial Reference-LAPACK form: no OpenMP or thread-private accumulators. +* + SUBROUTINE SLAED9_BR( K, NBR, D, QIN, LDQIN, QOUT, LDQOUT, + $ RHO, DLAMBDA, W, DELTA, WPROD, WSGN, + $ TAU, NEAR, INFO ) +* +* .. Scalar Arguments .. + INTEGER INFO, K, LDQIN, LDQOUT, NBR + REAL RHO +* .. +* .. Array Arguments .. + REAL D( * ), DELTA( * ), DLAMBDA( * ), + $ NEAR( 2, * ), QIN( LDQIN, * ), + $ QOUT( LDQOUT, * ), TAU( * ), W( * ), + $ WPROD( * ), WSGN( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) +* .. +* .. Local Scalars .. + INTEGER I, IORG, J, N1, N2, R + REAL ACC1, ACC2, DLT, ROWSUM, TEMP +* .. +* .. External Subroutines .. + EXTERNAL SLAED4, SLAED4_BR, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, SIGN, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 + IF( K.LT.0 ) THEN + INFO = -1 + ELSE IF( NBR.LT.0 ) THEN + INFO = -2 + ELSE IF( LDQIN.LT.MAX( 1, NBR ) ) THEN + INFO = -5 + ELSE IF( LDQOUT.LT.MAX( 1, NBR ) ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLAED9_BR', -INFO ) + RETURN + END IF +* + IF( K.EQ.0 ) + $ RETURN +* +* Root-only final merge: no parent boundary rows will be observed by +* a higher merge level, so updated W and secular vectors are not +* needed. +* + IF( NBR.EQ.0 ) THEN + DO 5 J = 1, K + CALL SLAED4( K, J, DLAMBDA, W, DELTA, RHO, D( J ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 170 + 5 CONTINUE + GO TO 170 + END IF +* +* For K = 1 or 2, SLAED4/SLAED5 returns normalized eigenvectors +* directly in DELTA. Keep this path identical to SLAED9's +* special-case copy to S. +* + IF( K.EQ.1 .OR. K.EQ.2 ) THEN + DO 30 J = 1, K + CALL SLAED4( K, J, DLAMBDA, W, DELTA, RHO, D( J ), + $ INFO ) + IF( INFO.NE.0 ) + $ GO TO 170 + IF( NBR.EQ.2 ) THEN + ACC1 = ZERO + ACC2 = ZERO + DO 10 I = 1, K + ACC1 = ACC1 + QIN( 1, I )*DELTA( I ) + ACC2 = ACC2 + QIN( 2, I )*DELTA( I ) + 10 CONTINUE + QOUT( 1, J ) = ACC1 + QOUT( 2, J ) = ACC2 + ELSE + DO 20 R = 1, NBR + ROWSUM = ZERO + DO 15 I = 1, K + ROWSUM = ROWSUM + QIN( R, I )*DELTA( I ) + 15 CONTINUE + QOUT( R, J ) = ROWSUM + 20 CONTINUE + END IF + 30 CONTINUE + GO TO 170 + END IF +* +* Preserve the update vector, solve all secular roots, then accumulate +* the product needed to reconstruct SLAED9's updated W. The +* high-relative-accuracy DELTA values produced by SLAED4_BR are +* consumed immediately rather than reconstructed later from the +* compact origin/tau representation. +* + DO 40 I = 1, K + WSGN( I ) = W( I ) + 40 CONTINUE +* + DO 41 I = 1, K + WPROD( I ) = ONE + 41 CONTINUE +* + DO 44 J = 1, K + CALL SLAED4_BR( K, J, DLAMBDA, WSGN, DELTA, RHO, D( J ), + $ IORG, TAU( J ), NEAR( 1, J ), + $ NEAR( 2, J ), INFO ) + IF( INFO.NE.0 ) + $ GO TO 170 +* + WPROD( J ) = WPROD( J )*DELTA( J ) + DO 42 I = 1, J - 1 + WPROD( I ) = WPROD( I )* + $ ( DELTA( I ) / + $ ( DLAMBDA( I )-DLAMBDA( J ) ) ) + 42 CONTINUE + DO 43 I = J + 1, K + WPROD( I ) = WPROD( I )* + $ ( DELTA( I ) / + $ ( DLAMBDA( I )-DLAMBDA( J ) ) ) + 43 CONTINUE + 44 CONTINUE +* + DO 80 I = 1, K + W( I ) = SIGN( SQRT( -WPROD( I ) ), WSGN( I ) ) + 80 CONTINUE +* +* Stream each secular eigenvector column through the selected rows. +* Reconstruct the DELTA column from the compact origin/tau state, +* using cached pole-adjacent values where cancellation is largest. +* + IF( NBR.EQ.2 ) THEN + DO 115 J = 1, K + IF( J.EQ.K ) THEN + IORG = K + N1 = K - 1 + N2 = K + ELSE IF( TAU( J ).GE.ZERO ) THEN + IORG = J + N1 = J + N2 = J + 1 + ELSE + IORG = J + 1 + N1 = J + N2 = J + 1 + END IF + TEMP = ZERO + ACC1 = ZERO + ACC2 = ZERO +* + DO 85 I = 1, N1 - 1 + DLT = W( I ) / + $ ( ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, I )*DLT + ACC2 = ACC2 + QIN( 2, I )*DLT + 85 CONTINUE +* + DLT = W( N1 ) / NEAR( 1, J ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, N1 )*DLT + ACC2 = ACC2 + QIN( 2, N1 )*DLT +* + DO 86 I = N1 + 1, N2 - 1 + DLT = W( I ) / + $ ( ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, I )*DLT + ACC2 = ACC2 + QIN( 2, I )*DLT + 86 CONTINUE +* + DLT = W( N2 ) / NEAR( 2, J ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, N2 )*DLT + ACC2 = ACC2 + QIN( 2, N2 )*DLT +* + DO 87 I = N2 + 1, K + DLT = W( I ) / + $ ( ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) ) + TEMP = TEMP + DLT*DLT + ACC1 = ACC1 + QIN( 1, I )*DLT + ACC2 = ACC2 + QIN( 2, I )*DLT + 87 CONTINUE +* + TEMP = SQRT( TEMP ) + QOUT( 1, J ) = ACC1 / TEMP + QOUT( 2, J ) = ACC2 / TEMP + 115 CONTINUE + GO TO 170 + END IF +* + DO 130 J = 1, K + IF( J.EQ.K ) THEN + IORG = K + ELSE IF( TAU( J ).GE.ZERO ) THEN + IORG = J + ELSE + IORG = J + 1 + END IF + TEMP = ZERO + DO 90 R = 1, NBR + QOUT( R, J ) = ZERO + 90 CONTINUE + DO 100 I = 1, K + DLT = ( DLAMBDA( I )-DLAMBDA( IORG ) ) - TAU( J ) + IF( J.EQ.K ) THEN + IF( I.EQ.K-1 ) + $ DLT = NEAR( 1, J ) + IF( I.EQ.K ) + $ DLT = NEAR( 2, J ) + ELSE + IF( I.EQ.J ) + $ DLT = NEAR( 1, J ) + IF( I.EQ.J+1 ) + $ DLT = NEAR( 2, J ) + END IF + DLT = W( I ) / DLT + TEMP = TEMP + DLT*DLT + DO 95 R = 1, NBR + QOUT( R, J ) = QOUT( R, J ) + QIN( R, I )*DLT + 95 CONTINUE + 100 CONTINUE + TEMP = SQRT( TEMP ) + DO 120 R = 1, NBR + QOUT( R, J ) = QOUT( R, J ) / TEMP + 120 CONTINUE + 130 CONTINUE +* + 170 CONTINUE + RETURN +* +* End of SLAED9_BR +* + END diff --git a/lapack-netlib/SRC/sstebr.f b/lapack-netlib/SRC/sstebr.f new file mode 100644 index 0000000000..be6604a7b4 --- /dev/null +++ b/lapack-netlib/SRC/sstebr.f @@ -0,0 +1,274 @@ +*> \brief \b SSTEBR +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE SSTEBR( N, D, E, WORK, LWORK, IWORK, LIWORK, INFO ) +* +* .. Scalar Arguments .. +* INTEGER INFO, LIWORK, LWORK, N +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL D( * ), E( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SSTEBR computes all eigenvalues of a symmetric tridiagonal matrix +*> using a values-only divide-and-conquer method that propagates only +*> the first and last boundary rows of each local eigenvector block +*> (STEBR = STymmetric tridiagonal Eigenvalues, Boundary-Row DC). +*> +*> Unlike SSTEDC with COMPZ = 'N', SSTEBR does not store full secular +*> eigenvector blocks across merge levels for later replay. It is a +*> distinct routine and does not replace SSTEDC or SSTERF. +*> +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix. N >= 0. +*> \endverbatim +*> +*> \param[in,out] D +*> \verbatim +*> D is REAL array, dimension (N) +*> On entry, the n diagonal elements of the tridiagonal matrix. +*> On exit, if INFO = 0, the eigenvalues in ascending order. +*> \endverbatim +*> +*> \param[in,out] E +*> \verbatim +*> E is REAL array, dimension (N-1) +*> On entry, the (n-1) subdiagonal elements of the tridiagonal +*> matrix. +*> On exit, E has been destroyed. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK. +*> If N <= SMLSIZ, where SMLSIZ is returned by ILAENV and is +*> typically about 25, LWORK must be at least 1. +*> If N > SMLSIZ, LWORK must be at least +*> MAX( 16*N, 2*N + SMLSIZ**2 + 2*SMLSIZ - 2 ). +*> +*> If LWORK = -1, then a workspace query is assumed; the +*> routine only calculates the optimal sizes of the WORK and +*> IWORK arrays, returns these values as the first entries of +*> the WORK and IWORK arrays, and no error message related to +*> LWORK or LIWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) +*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. +*> \endverbatim +*> +*> \param[in] LIWORK +*> \verbatim +*> LIWORK is INTEGER +*> The dimension of the array IWORK. +*> If N <= SMLSIZ, LIWORK must be at least 1. +*> If N > SMLSIZ, LIWORK must be at least 7*N. +*> +*> If LIWORK = -1, then a workspace query is assumed; the +*> routine only calculates the optimal sizes of the WORK and +*> IWORK arrays, returns these values as the first entries of +*> the WORK and IWORK arrays, and no error message related to +*> LWORK or LIWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: The algorithm failed to compute an eigenvalue while +*> working on a submatrix lying in rows and columns +*> INFO/(N+1) through MOD(INFO,N+1). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup auxOTHERcomputational +* +* ===================================================================== + SUBROUTINE SSTEBR( N, D, E, WORK, LWORK, IWORK, LIWORK, INFO ) +* +* -- LAPACK computational routine -- +* +* .. Scalar Arguments .. + INTEGER INFO, LIWORK, LWORK, N +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL D( * ), E( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) +* .. +* .. Local Scalars .. + LOGICAL LQUERY + INTEGER FINISH, LIWMIN, LWMIN, M, SMLSIZ, START + REAL EPS, ORGNRM, TINY +* .. +* .. External Functions .. + INTEGER ILAENV + REAL SLAMCH, SLANST, SROUNDUP_LWORK + EXTERNAL ILAENV, SLAMCH, SLANST, SROUNDUP_LWORK +* .. +* .. External Subroutines .. + EXTERNAL SLAED0_BR, SLASCL, SLASRT, SSTERF, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MOD, SQRT +* .. +* .. Executable Statements .. +* + INFO = 0 + LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) +* + IF( N.LT.0 ) THEN + INFO = -1 + END IF +* + IF( INFO.EQ.0 ) THEN + SMLSIZ = ILAENV( 9, 'SSTEBR', ' ', 0, 0, 0, 0 ) + IF( N.LE.SMLSIZ ) THEN + LWMIN = 1 + LIWMIN = 1 + ELSE +* +* BLO/BHI: 2*N; merge scratch: 14*N (see SLAED7_BR). +* Leaf solves need 2*N + SMLSIZ**2 + 2*SMLSIZ - 2: the +* SMLSIZ-by-SMLSIZ leaf eigenvector block plus the +* SSTEQR work array. +* + LWMIN = MAX( 16*N, 2*N + SMLSIZ*SMLSIZ + 2*SMLSIZ - 2 ) + LIWMIN = MAX( 1, 7*N ) + END IF + WORK( 1 ) = SROUNDUP_LWORK( LWMIN ) + IWORK( 1 ) = LIWMIN +* + IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN + INFO = -5 + ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SSTEBR', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + RETURN + END IF +* + IF( N.EQ.0 ) + $ RETURN + IF( N.EQ.1 ) + $ RETURN +* + IF( N.LE.SMLSIZ ) THEN + CALL SSTERF( N, D, E, INFO ) + GO TO 50 + END IF +* + ORGNRM = SLANST( 'M', N, D, E ) + IF( ORGNRM.LE.ZERO ) + $ GO TO 50 +* + EPS = SLAMCH( 'Epsilon' ) + START = 1 +* + 10 CONTINUE + IF( START.LE.N ) THEN + FINISH = START + 20 CONTINUE + IF( FINISH.LT.N ) THEN + TINY = EPS*SQRT( ABS( D( FINISH ) ) )* + $ SQRT( ABS( D( FINISH+1 ) ) ) + IF( ABS( E( FINISH ) ).GT.TINY ) THEN + FINISH = FINISH + 1 + GO TO 20 + END IF + END IF +* + M = FINISH - START + 1 + IF( M.EQ.1 ) THEN + START = FINISH + 1 + GO TO 10 + END IF +* + IF( M.GT.SMLSIZ ) THEN + ORGNRM = SLANST( 'M', M, D( START ), E( START ) ) + CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), + $ M, INFO ) + CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, + $ E( START ), M-1, INFO ) +* + CALL SLAED0_BR( M, D( START ), E( START ), WORK, IWORK, + $ INFO ) + IF( INFO.NE.0 ) THEN + INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) + + $ MOD( INFO, ( M+1 ) ) + START - 1 + GO TO 50 + END IF +* + CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), + $ M, INFO ) + ELSE + CALL SSTERF( M, D( START ), E( START ), INFO ) + IF( INFO.NE.0 ) THEN + INFO = START*( N+1 ) + FINISH + GO TO 50 + END IF + END IF +* + START = FINISH + 1 + GO TO 10 + END IF +* + CALL SLASRT( 'I', N, D, INFO ) +* + 50 CONTINUE + WORK( 1 ) = SROUNDUP_LWORK( LWMIN ) + IWORK( 1 ) = LIWMIN + RETURN +* +* End of SSTEBR +* + END diff --git a/lapack-netlib/TESTING/EIG/dchkst.f b/lapack-netlib/TESTING/EIG/dchkst.f index 6e02c84fbb..1e07cd142f 100644 --- a/lapack-netlib/TESTING/EIG/dchkst.f +++ b/lapack-netlib/TESTING/EIG/dchkst.f @@ -152,6 +152,9 @@ *> (26) | D1 - D2 | / ( |D1| ulp ) DSTEDC('V') and *> DSTEDC('N') *> +*> (38) | D1 - D5 | / ( |D1| ulp ) DSTEBR +*> (values-only boundary-row divide-and-conquer) +*> *> Test 27 is disabled at the moment because DSTEMR does not *> guarantee high relatvie accuracy. *> @@ -626,9 +629,9 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. Local Scalars .. LOGICAL BADNN, TRYRAC INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC, - $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, + $ JR, JSIZE, JSM, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, $ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS, - $ NMATS, NMAX, NSPLIT, NTEST, NTESTT + $ NMATS, NMAX, NSPLIT, NTEST, NTESTT, SMLSAV DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU @@ -636,7 +639,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), - $ KTYPE( MAXTYP ) + $ KTYPE( MAXTYP ), SMLSZS( 3 ) DOUBLE PRECISION DUMMA( 1 ) * .. * .. External Functions .. @@ -647,13 +650,15 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. External Subroutines .. EXTERNAL DCOPY, DLACPY, DLASET, DLASUM, DLATMR, DLATMS, $ DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD, DSTEBZ, - $ DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR, DSTERF, - $ DSTT21, DSTT22, DSYT21, DSYTRD, XERBLA + $ DSTEBR, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR, + $ DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, XERBLA, + $ XLAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT * .. * .. Data statements .. + DATA SMLSZS / 2, 8, 25 / DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8, $ 8, 8, 9, 9, 9, 9, 9, 10 / DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, @@ -1897,6 +1902,52 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ ULP*MAX( TEMP1, TEMP2 ) ) END IF 270 CONTINUE +* +* Call DSTEBR (values-only boundary-row DC), do test 38. +* Compare against D1 from DSTEQR('V'), same ratio as tests +* 12 and 26. Vary SMLSIZ to exercise merges even with the +* small matrix sizes in sep.in, and retain the worst ratio. +* + DO 265 JR = NTEST + 1, 37 + RESULT( JR ) = ZERO + 265 CONTINUE +* + NTEST = 38 + RESULT( 38 ) = ZERO + SMLSAV = ILAENV( 9, 'DSTEBR', ' ', 0, 0, 0, 0 ) + DO 267 JSM = 1, 3 + CALL DCOPY( N, SD, 1, D5, 1 ) + IF( N.GT.0 ) + $ CALL DCOPY( N-1, SE, 1, WORK, 1 ) +* + CALL XLAENV( 9, SMLSZS( JSM ) ) + CALL DSTEBR( N, D5, WORK, WORK( N+1 ), LWORK-N, + $ IWORK, LIWORK, IINFO ) + CALL XLAENV( 9, SMLSAV ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'DSTEBR', IINFO, N, + $ JTYPE, IOLDSD + INFO = ABS( IINFO ) + IF( IINFO.LT.0 ) THEN + RETURN + ELSE + RESULT( 38 ) = ULPINV + GO TO 280 + END IF + END IF +* + TEMP1 = ZERO + TEMP2 = ZERO +* + DO 266 J = 1, N + TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D5( J ) ) ) + TEMP2 = MAX( TEMP2, ABS( D1( J )-D5( J ) ) ) + 266 CONTINUE +* + RESULT( 38 ) = MAX( RESULT( 38 ), TEMP2 / + $ MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) ) + 267 CONTINUE +* 280 CONTINUE NTESTT = NTESTT + NTEST * diff --git a/lapack-netlib/TESTING/EIG/schkst.f b/lapack-netlib/TESTING/EIG/schkst.f index aecbdfe932..2b92f79004 100644 --- a/lapack-netlib/TESTING/EIG/schkst.f +++ b/lapack-netlib/TESTING/EIG/schkst.f @@ -152,6 +152,9 @@ *> (26) | D1 - D2 | / ( |D1| ulp ) SSTEDC('V') and *> SSTEDC('N') *> +*> (38) | D1 - D5 | / ( |D1| ulp ) SSTEBR +*> (values-only boundary-row divide-and-conquer) +*> *> Test 27 is disabled at the moment because SSTEMR does not *> guarantee high relatvie accuracy. *> @@ -626,9 +629,9 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. Local Scalars .. LOGICAL BADNN, TRYRAC INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC, - $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, + $ JR, JSIZE, JSM, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, $ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS, - $ NMATS, NMAX, NSPLIT, NTEST, NTESTT + $ NMATS, NMAX, NSPLIT, NTEST, NTESTT, SMLSAV REAL ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU @@ -636,7 +639,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), - $ KTYPE( MAXTYP ) + $ KTYPE( MAXTYP ), SMLSZS( 3 ) REAL DUMMA( 1 ) * .. * .. External Functions .. @@ -647,13 +650,15 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. External Subroutines .. EXTERNAL SCOPY, SLACPY, SLASET, SLASUM, SLATMR, SLATMS, $ SOPGTR, SORGTR, SPTEQR, SSPT21, SSPTRD, SSTEBZ, - $ SSTECH, SSTEDC, SSTEMR, SSTEIN, SSTEQR, SSTERF, - $ SSTT21, SSTT22, SSYT21, SSYTRD, XERBLA + $ SSTEBR, SSTECH, SSTEDC, SSTEMR, SSTEIN, SSTEQR, + $ SSTERF, SSTT21, SSTT22, SSYT21, SSYTRD, XERBLA, + $ XLAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN, REAL, SQRT * .. * .. Data statements .. + DATA SMLSZS / 2, 8, 25 / DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8, $ 8, 8, 9, 9, 9, 9, 9, 10 / DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, @@ -1897,6 +1902,52 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ ULP*MAX( TEMP1, TEMP2 ) ) END IF 270 CONTINUE +* +* Call SSTEBR (values-only boundary-row DC), do test 38. +* Compare against D1 from SSTEQR('V'), same ratio as tests +* 12 and 26. Vary SMLSIZ to exercise merges even with the +* small matrix sizes in sep.in, and retain the worst ratio. +* + DO 265 JR = NTEST + 1, 37 + RESULT( JR ) = ZERO + 265 CONTINUE +* + NTEST = 38 + RESULT( 38 ) = ZERO + SMLSAV = ILAENV( 9, 'SSTEBR', ' ', 0, 0, 0, 0 ) + DO 267 JSM = 1, 3 + CALL SCOPY( N, SD, 1, D5, 1 ) + IF( N.GT.0 ) + $ CALL SCOPY( N-1, SE, 1, WORK, 1 ) +* + CALL XLAENV( 9, SMLSZS( JSM ) ) + CALL SSTEBR( N, D5, WORK, WORK( N+1 ), LWORK-N, + $ IWORK, LIWORK, IINFO ) + CALL XLAENV( 9, SMLSAV ) + IF( IINFO.NE.0 ) THEN + WRITE( NOUNIT, FMT = 9999 )'SSTEBR', IINFO, N, + $ JTYPE, IOLDSD + INFO = ABS( IINFO ) + IF( IINFO.LT.0 ) THEN + RETURN + ELSE + RESULT( 38 ) = ULPINV + GO TO 280 + END IF + END IF +* + TEMP1 = ZERO + TEMP2 = ZERO +* + DO 266 J = 1, N + TEMP1 = MAX( TEMP1, ABS( D1( J ) ), ABS( D5( J ) ) ) + TEMP2 = MAX( TEMP2, ABS( D1( J )-D5( J ) ) ) + 266 CONTINUE +* + RESULT( 38 ) = MAX( RESULT( 38 ), TEMP2 / + $ MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) ) ) + 267 CONTINUE +* 280 CONTINUE NTESTT = NTESTT + NTEST *