Commit f006f8812a2cac3745df8d36b6eff980ed12a118

Authored by Aleksey
1 parent 23d187b7
Exists in master

DQDS

Showing 1 changed file with 121 additions and 0 deletions   Show diff stats
DQDS/dqds.c 0 → 100644
... ... @@ -0,0 +1,121 @@
  1 +#include <stdio.h>
  2 +#include <mkl.h>
  3 +#include "mpi.h"
  4 +
  5 +/* Auxiliary routines prototypes */
  6 +extern void print_matrix( char* desc, int m, int n, double* a, int lda );
  7 +
  8 +/* Parameters */
  9 +#define M 10
  10 +#define N M
  11 +#define LDA M
  12 +#define LDU M
  13 +#define LDVT N
  14 +
  15 +/* Main program */
  16 +int main(int argc,char* argv[]) {
  17 + /* Locals */
  18 + int size=100;
  19 + int myid,nprocs;
  20 + int i;
  21 +
  22 +
  23 + //std::cout<<argc<<" "<<argv[1]<<std::endl;
  24 +
  25 +
  26 + if (argc>1)
  27 + {
  28 + size = atoi(argv[1]);
  29 + }
  30 +
  31 +
  32 + MPI_Init(&argc,&argv);
  33 + MPI_Comm_rank(MPI_COMM_WORLD,&myid);
  34 + // int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info;
  35 + int m = size, n = size, lda = size, ldu = size, ldvt = size, info;
  36 + double wkopt;
  37 + /* Local arrays */
  38 + // double u[LDU*N];
  39 + // double s[N], u[LDU*N], vt[LDVT*N];
  40 + // double s[size], u[size*size], vt[size*size];
  41 + // double a[LDA*N] = {
  42 + // 8.79, 6.11, -9.15, 9.57, -3.49, 9.84,
  43 + // 9.93, 6.91, -793, 1.64, 4.02, 0.15,
  44 + // 9.83, 5.04, 4.86, 8.83, 9.80, -8.99,
  45 + // 9.83, 5.04, 4.86, 8.83, 9.80, -8.99,
  46 + // 9.83, 5.04, 4.86, 8.83, 9.80, -8.99,
  47 + // 9.83, 5.04, 4.86, 8.83, 9.80, -8.99,
  48 + // };
  49 + //
  50 + // double a[LDA*N];
  51 +
  52 + //double a[size*size];
  53 + // double* a = (double*) malloc(sizeof(double)*4);
  54 + // for (int i = 0; i < LDA*N; ++i)
  55 +
  56 + //for (int i = 0; i < size*size; i++)
  57 + //{
  58 + // a[i]=rand();
  59 + //}
  60 +
  61 + double diag[size],e_diag[size];
  62 + for ( i = 0; i<size;i++)
  63 + {
  64 + diag[i]=rand();
  65 + e_diag[i]=rand();
  66 + }
  67 + //e_diag[0]=0;
  68 +
  69 +
  70 + /* Executable statements */
  71 +
  72 + int ldc = 1;
  73 +
  74 + // info=LAPACKE_dgebrd(LAPACK_COL_MAJOR, m, n, a, lda, d, e, tauq, taup);
  75 + // info=LAPACKE_dgebrd(LAPACK_COL_MAJOR, m, n, a, lda, d, e, tauq, taup);
  76 + //printf( "DBDSQR Example Program Results\n" );
  77 + int ncvt = 0;
  78 + int nru = 0;
  79 + int ncc = 0;
  80 +
  81 + double starttime = MPI_Wtime();
  82 +
  83 + //info=LAPACKE_dorgbr( LAPACK_COL_MAJOR, 'Q', m, n, n, a, ldu, tauq);
  84 + //print_matrix( "Singular values", 1, n, diag, 1 );
  85 + info=LAPACKE_dbdsqr( LAPACK_COL_MAJOR, 'U', n, ncvt, nru, ncc, diag, e_diag,
  86 + NULL, ldvt, NULL, ldu, NULL, ldc);
  87 + //print_matrix( "Singular values", 1, n, diag, 1 );
  88 + if( info > 0 ) {
  89 + printf( "The algorithm computing SVD failed to converge.\n" );
  90 + exit( 1 );
  91 + }
  92 + /* Print singular values */
  93 + // print_matrix( "Singular values", 1, n, d, 1 );
  94 + /* Print left singular vectors */
  95 + // print_matrix( "Left singular vectors (stored columnwise)", m, n, u, ldu );
  96 + //delete [] d,e,;
  97 + if (!myid)
  98 + {
  99 + FILE *fp;
  100 + if ((fp=fopen("perf.log","a+"))==NULL)
  101 + {
  102 + printf("perf.log is unavalable");
  103 + }
  104 + fprintf(fp,"%d %d %f\n",1,size, MPI_Wtime() - starttime);
  105 + fclose(fp);
  106 + }
  107 +
  108 + print_matrix("Singular values are: ", 1, n, diag, 1);
  109 + MPI_Finalize();
  110 + exit(0);
  111 +} /* End of DBDSQR Example */
  112 +
  113 +/* Auxiliary routine: printing a matrix */
  114 +void print_matrix( char* desc, int m, int n, double* a, int lda ) {
  115 + int i, j;
  116 + printf( "\n %s\n", desc );
  117 + for( i = 0; i < m; i++ ) {
  118 + for( j = 0; j < n; j++ ) printf( " %12.8f", a[i+j*lda] );
  119 + printf( "\n" );
  120 + }
  121 +}
... ...