problem with LAPACKE_dstemr

A___Fiori — Mon, 11 Mar 2019 09:18:22 GMT

Hello! I am trying to use the function 'LAPACKE_dstemr' with a symmetric tridiagonal matrix 'T' but my code crashes. The matrix T (dimxdim) in general has the following form:

| 1.0 a 0 0 0 |
| a b a 0 0 |
T = | 0 a b a 0 | *tau
| 0 0 a b a |
| 0 0 0 a 1.0|

where a = -phi1, b=(1.0 + phi1^2)

Let for example that coefs=[phi1 tau]=[0.7 1] and dim=200. I want to compute all the eigenvalues and eigenvectors. I cant find an example with this function, please could you tell me what I am doing wrong. My mex-code is the following. Thank you very much.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include "mkl_vml.h"
#include "mex.h"
#include "matrix.h"
#include "mkl_vsl.h"
#include <string.h> /* is needed for the function 'strcmp' */

/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda );
extern double sprod(int dim, double *a);

/* main fucntion */
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
    int i,dim,info,n;
    int tryrac = true;
    double *coefs,phi1,tau,phi1square,rw2,*U,*w,*e,*d;
    
    /* make pointers to input data */
    dim = (int)mxGetScalar(prhs[0]);
    coefs = (double*)mxGetPr(prhs[1]);
    rw2 = (double)mxGetScalar(prhs[2]);
    
    /* phi1,tau */
    phi1 = coefs[0];
    phi1square = pow(phi1,2);
    tau = coefs[1];
    
    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix(1,dim,mxREAL);
    w = mxGetPr(plhs[0]);
    
    plhs[1] = mxCreateDoubleMatrix(dim,dim,mxREAL);
    U = mxGetPr(plhs[1]);
    
    /* Memory */
    d = (double*)mxMalloc(dim*sizeof(double));
    e = (double*)mxMalloc(dim*sizeof(double));
    
    /* e, d */
    e[0:dim:1] = -phi1*tau;
    d[0] = tau;
    d[dim-1] = tau;
    d[1:dim-2:1] = (1.0 + phi1square)*tau;
    
    print_matrix( "e", 1, dim, e, 1 );
    print_matrix( "d", 1, dim, d, 1 );
    
    /* Solve eigenproblem */
    n = dim;
    info = LAPACKE_dstemr(LAPACK_COL_MAJOR, 'V', 'A', dim,
            d, e, 0, 0, 0, 0, &n,  w, U, dim, dim, &n, &tryrac );
    
    mexPrintf( "n=%d\t dim=%d\n", n,dim );
    
    /* release memory */
    mxFree(e),mxFree(d);
    return ;
}

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda ) {
    MKL_INT i, j;
    mexPrintf( "\n %s\n", desc );
    for( i = 0; i < m; i++ ) {
        for( j = 0; j < n; j++ ) mexPrintf( " %6.2f", a[i+j*lda] );
        mexPrintf( "\n" );
    }
}

vl, vu

Gennady_F_Intel — Mon, 11 Mar 2019 13:17:44 GMT

vl, vu

If range = 'V', the lower and upper bounds of the interval to be searched for eigenvalues. Constraint: vl<vu.

topic vl, vu in Intel® oneAPI Math Kernel Library

problem with LAPACKE_dstemr

vl, vu