MKL LAPACK's dsptrd + dpteqr routines give (partly) wrong eigenvectors

MKL provides dsptrd and dpteqr routines to solve real symmetric positive-definite eigenvalue problems. dsptrd reduces a matrix to triangular form and dpteqr finds the eigenvalues and eigenvectors. But now I find that the eigenvectors given by the routines above are different from those given by dsyev, a routine that solves real symmetric (but not necessarily positive-definite) eigenvalue problems.

#include <iostream>
#include "mkl.h"

void RealSymPosiDefDiagonalization(double * matrix,int nrows,double ** eigenvalues,double ** eigenvectors){
    std::cout<<"MKL_LAPACK diagonalization for real symmetric positive-definite matrices with dsptrd and dpteqr... ";
    double * d=new double[nrows];*eigenvalues=d;
    double * e=new double[nrows-1];
    double * tau=new double[nrows-1];
    double * z=new double[nrows*nrows];*eigenvectors=z;
    MKL_INT info=LAPACKE_dsptrd(LAPACK_ROW_MAJOR,'L',nrows,matrix,d,e,tau);
    if (info==0){
        info=LAPACKE_dpteqr(LAPACK_ROW_MAJOR,'I',nrows,d,e,z,nrows);
        if (info==0){
            std::cout<<"Done"<<std::endl;
        }else{
            std::cout<<"Failed dpteqr_info="<<info<<std::endl;
        }
    }else{
        std::cout<<"Failed dsptrd_info="<<info<<std::endl;
    }
}

void RealSymDiagonalization(double * matrix,double * eigenvalues,int nrows){
    std::cout<<"MKL_LAPACK diagonalization for real symmetric matrices with dsyev ... ";
    MKL_INT info=LAPACKE_dsyev(LAPACK_ROW_MAJOR,'V','U',nrows,matrix,nrows,eigenvalues);
    if (info>0){
        std::cout<<"Failed info="<<info<<std::endl;
    }else{
        std::cout<<"Done"<<std::endl;
    }
}


void PrintMatrix(double * matrix,int nrows,int ncols,char storage){
    if (storage=='f'){
        for (int irow=0;irow<nrows;irow++){
            for (int icol=0;icol<ncols;icol++){
                std::cout<<" "<<matrix[irow*nrows+icol];
            }
            std::cout<<std::endl;
        }
    }else if (storage=='l'){
        for (int irow=0;irow<nrows;irow++){
            for (int icol=0;icol<ncols;icol++){
                if (irow>=icol){
                    std::cout<<" "<<matrix[irow*(irow+1)/2+icol];
                }else{
                    std::cout<<" "<<matrix[icol*(icol+1)/2+irow];
                }
            }
            std::cout<<std::endl;
        }
    }
}

int main(){
    double matrix[]={2  ,
                     1.5, 3  ,
                     1.5, 1.5, 5};
    std::cout<<"Target matrix:"<<std::endl;
    PrintMatrix(matrix,3,3,'l');
    double * eigenvalues;
    double * eigenvectors;
    RealSymPosiDefDiagonalization(matrix,3,&eigenvalues,&eigenvectors);
    std::cout<<"Eigenvalues: "<<std::endl;
    PrintMatrix(eigenvalues,1,3,'f');
    std::cout<<"Eigenvectors: "<<std::endl;
    PrintMatrix(eigenvectors,3,3,'f');


    double matrix1[]={2,1.5,1.5,
                      1.5,3,1.5,
                      1.5,1.5,5};
    RealSymDiagonalization(matrix1,eigenvalues,3);
    std::cout<<"Eigenvalues: "<<std::endl;
    PrintMatrix(eigenvalues,1,3,'f');
    std::cout<<"Eigenvectors: "<<std::endl;
    PrintMatrix(matrix1,3,3,'f'); // Now matrix1 is overwritten to be the eigenvector matrix.

    return 0;
}

The results are as follows.

Target matrix:
 2 1.5 1.5
 1.5 3 1.5
 1.5 1.5 5
MKL_LAPACK diagonalization for real symmetric positive-definite matrices with dsptrd and dpteqr... Done
Eigenvalues: 
 6.69646 2.42659 0.876946
Eigenvectors: 
 0.402264 -0.342101 -0.849206
 -0.890587 0.0687953 -0.44958
 0.212223 0.937142 -0.276997
MKL_LAPACK diagonalization for real symmetric matrices with dsyev ... Done
Eigenvalues: 
 0.876946 2.42659 6.69646
Eigenvectors: 
 0.849206 0.342101 0.402264
 -0.513767 0.711305 0.479675
 -0.122035 -0.614014 0.779804

The two methods give identical eigenvalues and first rows of the eigenvectors, but their second and third rows of the eigenvectors are different. As far as I know, dsyev has given the true answer. Very strange. Does anyone have any idea? Thanks for any help!