4. Matrixrang - SVD (GSL) - heap allocation

Matrixrang - SVD (GSL) - heap allocation

  • 6

    Hallo,

    ich versuche mich gerade an der Berechnung des Ranges einer Matrix und verwende dazu GSL.
    Nun erhalte ich bei folgendem Programm

    #include "utils.h"

/* Computes the rank of a matrix using Singular Value Decomposition */ 

long rank(long n_rows, long n_cols, double **M)
{
  long row, col, i, col_row_min, rank = 0, transp = 0; 

  double **U, *S/*, *X*/, *V, *WORK;

  if ((n_rows <=0) || (n_cols <= 0))
    return -1;		

  col_row_min = (n_cols < n_rows ? n_cols : n_rows);

  U = (double **) malloc (n_rows * n_cols * sizeof (double));
  if (!U)
  {
    perror ("__FILE__,__LINE__");

    return -1;
  }

  S = (double*) malloc (col_row_min * sizeof(double));
  if (!S)
  {
    perror ("__FILE__,__LINE__");
    free ((void*)U);

    return -1;
  }

  V = (double*) malloc (col_row_min * col_row_min * sizeof(double));
  if (!V)
  {
    perror ("__FILE__,__LINE__");
    free ((void*)U);
    free ((void*)S);

    return -1;
  }

/*  X = (double*) malloc (col_row_min * col_row_min * sizeof(double));
  if (!X)
  {
    perror ("__FILE__,__LINE__");
    free ((void*)U);
    free ((void*)V);
    free ((void*)S);

    return -1;
  }
*/
  WORK = (double*) malloc (col_row_min  * sizeof(double));
  if (!WORK)
  {
    perror ("__FILE__,__LINE__");
    free ((void*)U);
    free ((void*)V);
    free ((void*)S);
//  free ((void*)X);

    return -1;
  }

  if (n_rows<n_cols) 
  {
    for (row = 0; row < n_rows; row++)
      for (col = 0; col < n_cols; col++)
        U[row][col] = M[row][col];
  }
  else
  {
    transp = 1;
    for (row = 0; row < n_rows; row++)
      for (col = 0; col < n_cols; col++)
        U[col][row] = M[row][col];
  }

  //gsl_linalg_SV_decomp_mod (&U, &X, &V, &S, &WORK);
  gsl_linalg_SV_decomp(&U, &V, &S, &WORK);

  if (_ABS(S[0])>EPSILON)
    rank = 1;
  else
    rank = 0;
  for (i=1; i< col_row_min;i++)
    if (_ABS(S[i])>EPSILON)
      rank ++;

  free((void*)U);
  free((void*)V);
  free((void*)S);
//free((void*)X);
  free((void*)WORK);

  return rank;
}

    mit valgrind die Meldung

    ==7464== Conditional jump or move depends on uninitialised value(s)
==7464==    at 0x4EDB104: gsl_linalg_SV_decomp (svd.c:63)
==7464==    by 0x4014BF: rank (utils.c:83)
[...]
==7464==  Uninitialised value was created by a heap allocation
==7464==    at 0x4C2815C: malloc (vg_replace_malloc.c:236)
[...]
==7464== 
==7464== Conditional jump or move depends on uninitialised value(s)
==7464==    at 0x4EDB194: gsl_linalg_SV_decomp (svd.c:68)
==7464==    by 0x4014BF: rank (utils.c:83)
[...]
==7464==  Uninitialised value was created by a heap allocation
==7464==    at 0x4C2815C: malloc (vg_replace_malloc.c:236)
[...]

    Sieht jemand einen Fehler?

    Gruß,
    6H05T

    0
  • Mod

    6H05T schrieb:

    Sieht jemand einen Fehler?

    Wenn du dein Debugging an das Forum auslagern willst, dann liefer auch einen debugbaren (=compilierbaren) Code.

    0
  • B

    Laut Doku benötigt gsl_linalg_SV_decomp Zeiger auf gsl_matrix. Du stopfst einfach double-Zeiger rein. Dein Compiler sollte davor eigentlich warnen.

    0
  • 6

    SeppJ schrieb:

    6H05T schrieb:

    Sieht jemand einen Fehler?

    Wenn du dein Debugging an das Forum auslagern willst, dann liefer auch einen debugbaren (=compilierbaren) Code.

    Nein, das war eigentlich nicht meine Idee.

    Bashar schrieb:

    Laut Doku benötigt gsl_linalg_SV_decomp Zeiger auf gsl_matrix. Du stopfst einfach double-Zeiger rein. Dein Compiler sollte davor eigentlich warnen.

    Ja, danke 👍 . Das hatte ich übersehen.

    Mein nächstes Problem ist nun, dass ich nicht genau weiß, wie ich die Matrix korrekt an gsl_linalg_SV_decomp übergebe (gsl_matrix verwenden ??).

    0
  • B

    6H05T schrieb:

    Mein nächstes Problem ist nun, dass ich nicht genau weiß, wie ich die Matrix korrekt an gsl_linalg_SV_decomp übergebe (gsl_matrix verwenden ??).

    Ja natürlich gsl_matrix verwenden. Lies die Dokumentation zur GSL. Da gibt es sogar Beispiele.

    0
  • 6

    Ja, natürlich. Dumme Frage.
    Muss mal eben sehen, wie ich die einzelnen Einträge von M an eine gsl_matrix übergebe.

    0
  • 6

    Super, wie einfach das ist.
    Alles sofort aus der Dok zu erkennen.

    0
  • 6

    Muss ich für gsl_matrix und gsl_vector noch zusätzlich etwas einbinden.
    Beim kompilieren erkennt er diese nämlich nicht.
    gsl_matrix.h und gsl_vector.h?
    Wären die korrekt?

    long row, col, i, col_row_min, col_row_max, rank = 0, transp = 0; 

  gsl_matrix *U, *V;
  gsl_vector *S, *WORK;

  if ((n_rows <=0) || (n_cols <= 0))
    return -1;		

  col_row_min = (n_cols < n_rows ? n_cols : n_rows);

  S = gsl_vector_alloc (col_row_min);
  WORK = gsl_vector_alloc (col_row_min);
  V = gsl_matrix_alloc (n_rows,n_cols);
  U = gsl_matrix_alloc (n_rows,n_cols);
/*  Fehlerüberprüfung kommt noch!!! */

  if (n_rows<n_cols) 
  {
    V = gsl_matrix_alloc (n_rows,n_cols);
    U = gsl_matrix_alloc (n_rows,n_cols);
    for (row = 0; row < n_rows; row++)
      for (col = 0; col < n_cols; col++)
        gsl_matrix_set(U,row,col,M[row][col]);
  }
  else
  {
    V = gsl_matrix_alloc (n_cols,n_rows);
    U = gsl_matrix_alloc (n_cols,n_rows);
    transp = 1;
    for (row = 0; row < n_rows; row++)
      for (col = 0; col < n_cols; col++)
        gsl_matrix_set(U,row,col,M[col][row]);
  }

  gsl_linalg_SV_decomp(U, V, S, WORK);

  if (_ABS(S[0])>EPSILON)
    rank = 1;
  else
    rank = 0;
  for (i=1; i< col_row_min;i++)
    if (_ABS(S[i])>EPSILON)
      rank ++;

  gsl_vector_free(WORK);
  gsl_vector_free(S);
  gsl_matrix_free(U);
  gsl_matrix_free(V);

  return rank;
}
    0
  • ?

    Aus 8.4 der Doku:

    GSL Docu schrieb:

    The functions for allocating and accessing matrices are defined in gsl_matrix.h

    2.2 Compiling and Linking:

    GSL Docu schrieb:

    The library header files are installed in their own gsl directory. You should write any preprocessor include statements with a gsl/ directory prefix thus,

    #include <gsl/gsl_math.h>

    0
Anmelden zum Antworten