Justine Devin / gammalib

Commit 983adf00a29da47c81edf5d69cfeda45cb43eb36

Authored by knodlseder
1 parent bf27cf6f

correct bug in cs_transpose as spotted by ISDC guys

Inline Side-by-side
Showing 1 changed file with 129 additions and 103 deletions
src/linalg/GSparseMatrix.cpp
  View file @983adf0
... ... @@ -4,19 +4,32 @@
4 4 * copyright (C) 2006-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7   - * This program is free software; you can redistribute it and/or modify *
8   - * it under the terms of the GNU General Public License as published by *
9   - * the Free Software Foundation; either version 2 of the License, or *
10   - * (at your option) any later version. *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
11 19 * *
12 20 * ----------------------------------------------------------------------- *
13 21 * This class implements the compressed sparse column format. The *
14 22 * following arrays are allocated: *
15   - * m_data : holds all 'elements' non-zero values, in column order (Ax) *
16   - * m_colstart : holds the index of the first element of each column (Ap) *
17   - * m_rowinx : holds the row indices for all elements (Ai) *
  23 + * *
  24 + * m_rows (n) Number of rows *
  25 + * m_cols (m) Number of columns *
  26 + * m_data (Ax) Holds all 'elements' non-zero values, in column order *
  27 + * m_colstart (Ap) Holds the index of the first element of each column *
  28 + * m_rowinx (Ai) Holds the row indices for all elements *
  29 + * *
18 30 * Column 'col' covers therefore [m_colstart[j], ..., m_colstart[j+1]-1] *
19 31 * Note that 'm_colstart' has m_cols+1 elements. *
  32 + * *
20 33 ***************************************************************************/
21 34 /**
22 35 * @file GSparseMatrix.cpp
... ... @@ -2829,17 +2842,17 @@ GSparseMatrix abs(const GSparseMatrix&amp; m)
2829 2842 ***************************************************************************/
2830 2843 GSparseMatrix cs_symperm(const GSparseMatrix& m, const int* pinv)
2831 2844 {
2832   - // Declare loop variables
2833   - int i, j, p, q, i2, j2;
  2845 + // Declare loop variables
  2846 + int i, j, p, q, i2, j2;
2834 2847  
2835   - // Assign matrix attributes
2836   - int n = m.m_cols;
2837   - int* Ap = m.m_colstart;
2838   - int* Ai = m.m_rowinx;
2839   - double* Ax = m.m_data;
  2848 + // Assign matrix attributes
  2849 + int n = m.m_cols;
  2850 + int* Ap = m.m_colstart;
  2851 + int* Ai = m.m_rowinx;
  2852 + double* Ax = m.m_data;
2840 2853  
2841   - // Allocate result matrix
2842   - GSparseMatrix C(n, n, Ap[n]);
  2854 + // Allocate result matrix
  2855 + GSparseMatrix C(n, n, Ap[n]);
2843 2856  
2844 2857 // Allocate and initialise workspace
2845 2858 int wrk_size = n;
... ... @@ -2847,66 +2860,69 @@ GSparseMatrix cs_symperm(const GSparseMatrix&amp; m, const int* pinv)
2847 2860 for (i = 0; i < wrk_size; ++i)
2848 2861 wrk_int[i] = 0;
2849 2862  
2850   - // Assign result matrix attributes
2851   - int* Cp = C.m_colstart;
2852   - int* Ci = C.m_rowinx;
2853   - double* Cx = C.m_data;
  2863 + // Assign result matrix attributes
  2864 + int* Cp = C.m_colstart;
  2865 + int* Ci = C.m_rowinx;
  2866 + double* Cx = C.m_data;
2854 2867  
2855   - // Count entries in each column of C
2856   - for (j = 0; j < n; j++) {
  2868 + // Count entries in each column of C
  2869 + for (j = 0; j < n; j++) {
2857 2870  
2858   - // Column j of A is column j2 of C
2859   - j2 = pinv ? pinv[j] : j;
  2871 + // Column j of A is column j2 of C
  2872 + j2 = pinv ? pinv[j] : j;
2860 2873  
2861   - // Loop over entries in column j
2862   - for (p = Ap[j]; p < Ap[j+1]; p++) {
2863   - i = Ai [p];
2864   - if (i > j) continue; // skip lower triangular part of A
2865   - i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
2866   - wrk_int[G_MAX(i2, j2)]++; // column count of C
  2874 + // Loop over entries in column j
  2875 + for (p = Ap[j]; p < Ap[j+1]; p++) {
  2876 + i = Ai [p];
  2877 + if (i > j) continue; // skip lower triangular part of A
  2878 + i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
  2879 + wrk_int[G_MAX(i2, j2)]++; // column count of C
  2880 + }
2867 2881 }
2868   - }
2869 2882  
2870   - // Compute column pointers of C
2871   - cs_cumsum(Cp, wrk_int, n);
  2883 + // Compute column pointers of C
  2884 + cs_cumsum(Cp, wrk_int, n);
2872 2885  
2873   - // Loop over all columns of A
2874   - for (j = 0 ; j < n ; j++) {
  2886 + // Loop over all columns of A
  2887 + for (j = 0 ; j < n ; j++) {
2875 2888  
2876   - // Column j of A is column j2 of C
2877   - j2 = pinv ? pinv[j] : j;
  2889 + // Column j of A is column j2 of C
  2890 + j2 = pinv ? pinv[j] : j;
2878 2891  
2879   - // Loop over entries in column j
2880   - for (p = Ap[j]; p < Ap[j+1]; p++) {
2881   - i = Ai [p] ;
2882   - if (i > j) continue; // skip lower triangular part of A
2883   - i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
2884   - Ci[q = wrk_int[G_MAX(i2,j2)]++] = G_MIN(i2,j2);
2885   - if (Cx) Cx[q] = Ax[p];
  2892 + // Loop over entries in column j
  2893 + for (p = Ap[j]; p < Ap[j+1]; p++) {
  2894 + i = Ai [p] ;
  2895 + if (i > j) continue; // skip lower triangular part of A
  2896 + i2 = pinv ? pinv[i] : i; // row i of A is row i2 of C
  2897 + Ci[q = wrk_int[G_MAX(i2,j2)]++] = G_MIN(i2,j2);
  2898 + if (Cx) Cx[q] = Ax[p];
  2899 + }
2886 2900 }
2887   - }
2888 2901  
2889   - // Free workspace
2890   - delete [] wrk_int;
  2902 + // Free workspace
  2903 + delete [] wrk_int;
2891 2904  
2892   - // Rectify the number of elements in matrix C
2893   - C.free_elements(Cp[n], (C.m_elements-Cp[n]));
  2905 + // Rectify the number of elements in matrix C
  2906 + C.free_elements(Cp[n], (C.m_elements-Cp[n]));
2894 2907  
2895   - // Return result
2896   - return C;
  2908 + // Return result
  2909 + return C;
2897 2910  
2898 2911 }
2899 2912  
2900 2913  
2901 2914 /***************************************************************************
2902   - * cs_transpose *
2903   - * ----------------------------------------------------------------------- *
2904   - * Transpose matrix. The flag 'values' allows to avoid copying the actual *
2905   - * data values. This allows to perform a logical matrix transposition, as *
2906   - * needed by the symbolic matrix analysis class. *
2907   - * ----------------------------------------------------------------------- *
2908   - * NOTE: This routine DOES NOT support pending elements (they have to be *
2909   - * filled before. *
  2915 + * @brief Compute transpose matrix
  2916 + *
  2917 + * @param[in] m Matrix.
  2918 + * @param[in] value Flag that signals if values should be copied.
  2919 + *
  2920 + * The flag 'values' allows to avoid copying the actual data values. This
  2921 + * allows to perform a logical matrix transposition, as needed by the
  2922 + * symbolic matrix analysis class.
  2923 + *
  2924 + * Note that this method does not support pending elements (they have to be
  2925 + * filled before).
2910 2926 ***************************************************************************/
2911 2927 GSparseMatrix cs_transpose(const GSparseMatrix& m, int values)
2912 2928 {
... ... @@ -2916,45 +2932,55 @@ GSparseMatrix cs_transpose(const GSparseMatrix&amp; m, int values)
2916 2932 // Allocate and initialise workspace
2917 2933 int wrk_size = m.m_rows;
2918 2934 int* wrk_int = new int[wrk_size];
2919   - for (int i = 0; i < wrk_size; i++)
  2935 + for (int i = 0; i < wrk_size; ++i) {
2920 2936 wrk_int[i] = 0;
  2937 + }
2921 2938  
2922   - // Setup the number of non-zero elements in each row
2923   - for (int p = 0; p < m.m_elements; p++)
2924   - wrk_int[m.m_rowinx[p]]++;
2925   -
2926   - // Set row pointers. To use a GSparseSymbolic function we have to
2927   - // allocate and object (but this does not take memory)
2928   - cs_cumsum(result.m_colstart, wrk_int, m.m_rows);
2929   -
2930   - // Case A: Normal transponse, including assignment of values
2931   - if (values) {
2932   - for (int col = 0; col < m.m_cols; col++) {
2933   - for (int p = m.m_colstart[col] ; p < m.m_colstart[col+1] ; p++) {
2934   - int i = wrk_int[m.m_rowinx[p]]++;
2935   - result.m_rowinx[i] = col;
2936   - result.m_data[i] = m.m_data[p] ;
2937   - }
  2939 + // Setup the number of non-zero elements in each row
  2940 + // for (p = 0 ; p < Ap [n] ; p++) w [Ai [p]]++ ;
  2941 + // Ap[n] = m.m_colstart[col]
  2942 + // n = m.m_cols
  2943 + // Ai[p] = m.m_rowinx[p]
  2944 + for (int p = 0; p < m.m_colstart[m.m_cols]; ++p) {
  2945 + wrk_int[m.m_rowinx[p]]++;
2938 2946 }
2939   - }
2940 2947  
2941   - // Case B: Logical transponse, no assignment of values is performed
2942   - else {
2943   - for (int col = 0; col < m.m_cols; col++) {
2944   - for (int p = m.m_colstart[col] ; p < m.m_colstart[col+1] ; p++)
2945   - result.m_rowinx[wrk_int[m.m_rowinx[p]]++] = col;
  2948 + // Set row pointers. To use a GSparseSymbolic function we have to
  2949 + // allocate and object (but this does not take memory)
  2950 + cs_cumsum(result.m_colstart, wrk_int, m.m_rows);
  2951 +
  2952 + // Case A: Normal transponse, including assignment of values
  2953 + if (values) {
  2954 + for (int col = 0; col < m.m_cols; ++col) {
  2955 + for (int p = m.m_colstart[col]; p < m.m_colstart[col+1] ; ++p) {
  2956 + int i = wrk_int[m.m_rowinx[p]]++;
  2957 + result.m_rowinx[i] = col;
  2958 + result.m_data[i] = m.m_data[p] ;
  2959 + }
  2960 + }
2946 2961 }
2947   - }
2948 2962  
2949   - // Return transponse matrix
2950   - return result;
  2963 + // Case B: Logical transponse, no assignment of values is performed
  2964 + else {
  2965 + for (int col = 0; col < m.m_cols; ++col) {
  2966 + for (int p = m.m_colstart[col]; p < m.m_colstart[col+1] ; ++p) {
  2967 + result.m_rowinx[wrk_int[m.m_rowinx[p]]++] = col;
  2968 + }
  2969 + }
  2970 + }
  2971 +
  2972 + // Free workspace
  2973 + delete [] wrk_int;
  2974 +
  2975 + // Return transponse matrix
  2976 + return result;
2951 2977  
2952 2978 }
2953 2979  
2954 2980  
2955 2981 /*==========================================================================
2956 2982 = =
2957   - = Other functions used by GSparseMatrix =
  2983 + = Other functions =
2958 2984 = =
2959 2985 ==========================================================================*/
2960 2986  
... ... @@ -2970,22 +2996,22 @@ GSparseMatrix cs_transpose(const GSparseMatrix&amp; m, int values)
2970 2996 ***************************************************************************/
2971 2997 double cs_cumsum(int* p, int* c, int n)
2972 2998 {
2973   - // Signal error if one of the input pointers is NULL
2974   - if (!p || !c) return (-1);
2975   -
2976   - // Initialise sums (integer and double)
2977   - int nz = 0;
2978   - double nz2 = 0.0;
2979   -
2980   - // Evaluate p[0..n] = cumulative sum of c[0..n-1]
2981   - for (int i = 0; i < n; ++i) {
2982   - p[i] = nz ;
2983   - nz += c[i];
2984   - nz2 += c[i]; // also in double to avoid int overflow
2985   - c[i] = p[i]; // also copy p[0..n-1] back into c[0..n-1]
2986   - }
2987   - p[n] = nz ;
  2999 + // Signal error if one of the input pointers is NULL
  3000 + if (!p || !c) return (-1);
  3001 +
  3002 + // Initialise sums (integer and double)
  3003 + int nz = 0;
  3004 + double nz2 = 0.0;
  3005 +
  3006 + // Evaluate p[0..n] = cumulative sum of c[0..n-1]
  3007 + for (int i = 0; i < n; ++i) {
  3008 + p[i] = nz ;
  3009 + nz += c[i];
  3010 + nz2 += c[i]; // also in double to avoid int overflow
  3011 + c[i] = p[i]; // also copy p[0..n-1] back into c[0..n-1]
  3012 + }
  3013 + p[n] = nz ;
2988 3014  
2989   - // Return cumulative sum of c[0..n-1]
2990   - return nz2;
  3015 + // Return cumulative sum of c[0..n-1]
  3016 + return nz2;
2991 3017 }
... ...