Singular value decomposition to solve the null space problem

scd.cpp

#include "mkl.h"
#include <math.h>
#include <float.h>
#include <ctime>
#include <stdio.h>

/* Auxiliary routine: printing a matrix */
void print_matrix(double * matrix, int nrow, int ncol) {
  for (int i = 0; i < nrow; i++) {
    for (int j = 0; j < ncol; j++) {
      //printf("%3.1f\t",matrix[i * ncol + j]);
      printf("%.*e\t",FLT_DIG, matrix[i * ncol + j]);
    }
    printf("\n");
  }
}
/* Auxiliary routine: printing a matrix */
void octave_matrix(double * matrix, int nrow, int ncol) {
    printf("[");
  for (int i = 0; i < nrow; i++) {
    for (int j = 0; j < ncol; j++) {
      printf("%.*e ",FLT_DIG, matrix[i * ncol + j]);
    }
    printf(";");
  }
    printf("]\n");
}
int main () {
    int debug=0;
    int size = 64;
    for (int m = 3; m < 600; m ++) {
    double * A = (double *)mkl_calloc(m*m,sizeof(double),size);
    double * Acopy = (double *)mkl_calloc(m*m,sizeof(double),size);
    double * b = (double *) mkl_calloc(m,sizeof(double),size);
    double * z = (double *) mkl_calloc(m,sizeof(double),size);
    double * u = (double*)mkl_calloc(m*m,sizeof(double),size);
    double * s = (double*)mkl_calloc(m,sizeof(double),size);
    double * superb = (double*)mkl_calloc(m*m,sizeof(double),size);
    srand(unsigned(time(0)));
    for (int i = 0; i < m; i++) {
        double rowsum = 0;
        for (int j = 0; j < m; j++) {
            A[i*m+j] = (double) rand() / (double) (RAND_MAX-1);
            rowsum += A[i*m+j];
        }
        A[i*m+i]-=rowsum;
    }
    if (debug) {
        printf("A=\n");
        print_matrix(A,m,m);
    }

    mkl_domatcopy('R','N',m,m,1,A,m,Acopy,m);
    int matrix_order = LAPACK_ROW_MAJOR;
    char jobu = 'A';
    char jobvt = 'N';

    int info;
    
    info = LAPACKE_dgesvd(matrix_order,
            jobu, jobvt,
            m,m,
            Acopy,m,
            s, u, m,
            NULL, m,
            superb);

    for (int i = 0; i < m; i ++ ) {
        b[i]=u[i*m+(m-1)];
    }

    if (debug) {
        printf("$?=%d\nA=\n",info);
        print_matrix(A,m,m);
        octave_matrix(A,m,m);
        printf("U=\n");
        print_matrix(u,m,m);
        printf("b=\n");
        octave_matrix(b,m,1);
    }
    cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,m,1,m,1,A,m,b,1,1,z,1);

    if (debug) {
        printf("Z=\n");
        print_matrix(z,m,1);
    }
    double sum = 0;
    for (int i = 0; i < m; i ++) {
        sum+=fabs(z[i]);
    }
    printf("%d: %3.7f\n",m,sum);

    /*mkl_free(A);
    mkl_free(Acopy);
    //mkl_free(b);
    mkl_free(u);
    mkl_free(s);
    mkl_free(superb);*/
    }

    return 0;
}