21378561283

Committed 26 Jan 2026 11:50PM UTC coverage: 84.938% (-4.1%) from 88.988%

Build # 21378561283

Build Type

Pull #33

github

Committed by

web-flow

Commit Message

Merge f819d76ac into 8137ce2f7

Pull Request Pull Request #33: Windows build

Coverage Stats

1191 of 1400 branches covered (85.07%)

Branch coverage included in aggregate %.

149 of 158 new or added lines in 19 files covered. (94.3%)

171 existing lines in 11 files now uncovered.

3687 of 4343 relevant lines covered (84.9%)

130.32 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

82.55

/src/core/Matrix.c

/*
 * Copyright 2025 Daniel Cederberg
 *
 * This file is part of the PSLP project (LP Presolver).
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "Matrix.h"
#include "Debugger.h"
#include "Memory_wrapper.h"
#include "Numerics.h"
#include "RowColViews.h"
#include "glbopts.h"
#include "stdlib.h"
#include "string.h"

Matrix *matrix_new(const double *Ax, const int *Ai, const int *Ap, int n_rows,
                   int n_cols, int nnz)
{
    DEBUG(ASSERT_NO_ZEROS_D(Ax, ((size_t) (nnz))));
    Matrix *A = matrix_alloc(n_rows, n_cols, nnz);
    RETURN_PTR_IF_NULL(A, NULL);
    int offset, i, row_size, row_alloc;

    offset = 0;
    for (i = 0; i < n_rows; ++i)
    {
        A->p[i].start = Ap[i] + offset;
        memcpy(A->x + A->p[i].start, Ax + Ap[i],
               ((size_t) (Ap[i + 1] - Ap[i])) * sizeof(double));

        memcpy(A->i + A->p[i].start, Ai + Ap[i],
               ((size_t) (Ap[i + 1] - Ap[i])) * sizeof(int));

        A->p[i].end = Ap[i + 1] + offset;
        row_size = A->p[i].end - A->p[i].start;
        row_alloc = calc_memory_row(row_size, EXTRA_ROW_SPACE, EXTRA_MEMORY_RATIO);
        offset += row_alloc - row_size;
    }

    A->p[n_rows].start = Ap[n_rows] + offset;
    A->p[n_rows].end = A->p[n_rows].start;

    return A;
}

// needed for the transpose function
Matrix *matrix_alloc(int n_rows, int n_cols, int nnz)
{
    Matrix *A = (Matrix *) ps_malloc(1, sizeof(Matrix));
    RETURN_PTR_IF_NULL(A, NULL);

    A->m = n_rows;
    A->n = n_cols;
    A->nnz = nnz;
    A->n_alloc = calc_memory(nnz, n_rows, EXTRA_ROW_SPACE, EXTRA_MEMORY_RATIO);

#ifdef TESTING
    A->i = (int *) ps_calloc(((size_t) (A->n_alloc)), sizeof(int));
    A->p = (RowRange *) ps_calloc(((size_t) (n_rows + 1)), sizeof(RowRange));
    A->x = (double *) ps_calloc(((size_t) (A->n_alloc)), sizeof(double));
#else
    A->i = (int *) ps_malloc((size_t) (A->n_alloc), sizeof(int));
    A->p = (RowRange *) ps_malloc((size_t) (n_rows + 1), sizeof(RowRange));
    A->x = (double *) ps_malloc((size_t) (A->n_alloc), sizeof(double));
#endif

    if (!A->i || !A->p || !A->x)
    {
        free_matrix(A);
        return NULL;
    }

    return A;
}

Matrix *matrix_new_no_extra_space(const double *Ax, const int *Ai, const int *Ap,
                                  int n_rows, int n_cols, int nnz)
{
    Matrix *A = (Matrix *) ps_malloc(1, sizeof(Matrix));
    RETURN_PTR_IF_NULL(A, NULL);

    A->m = n_rows;
    A->n = n_cols;
    A->nnz = nnz;
    A->n_alloc = nnz;
    A->i = (int *) ps_malloc((size_t) A->n_alloc, sizeof(int));
    A->p = (RowRange *) ps_malloc((size_t) (n_rows + 1), sizeof(RowRange));
    A->x = (double *) ps_malloc((size_t) A->n_alloc, sizeof(double));

    if (!A->i || !A->p || !A->x)
    {
        free_matrix(A);
        return NULL;
    }

    memcpy(A->x, Ax, ((size_t) (nnz)) * sizeof(double));
    memcpy(A->i, Ai, ((size_t) (nnz)) * sizeof(int));

    for (int i = 0; i <= n_rows; ++i)
    {
        A->p[i].start = Ap[i];
        A->p[i].end = Ap[i + 1];
    }

    return A;
}

Matrix *transpose(const Matrix *A, int *work_n_cols)
{
    Matrix *AT = matrix_alloc(A->n, A->m, A->nnz);
    RETURN_PTR_IF_NULL(AT, NULL);
    int i, j, start;
    int *count = work_n_cols;
    memset(count, 0, ((size_t) (A->n)) * sizeof(int));

    // -------------------------------------------------------------------
    //  compute nnz in each column of A
    // -------------------------------------------------------------------
    for (i = 0; i < A->m; ++i)
    {
        for (j = A->p[i].start; j < A->p[i].end; ++j)
        {
            count[A->i[j]]++;
        }
    }
    // ------------------------------------------------------------------
    //  compute row pointers, taking the extra space into account
    // ------------------------------------------------------------------
    AT->p[0].start = 0;
    for (i = 0; i < A->n; ++i)
    {
        start = AT->p[i].start;
        AT->p[i].end = start + count[i];
        AT->p[i + 1].start =
            start + calc_memory_row(count[i], EXTRA_ROW_SPACE, EXTRA_MEMORY_RATIO);
        count[i] = start;
    }

    AT->p[A->n].start = AT->n_alloc;
    AT->p[A->n].end = AT->n_alloc;

    // ------------------------------------------------------------------
    //  fill transposed matrix (this is a bottleneck)
    // ------------------------------------------------------------------
    for (i = 0; i < A->m; ++i)
    {
        for (j = A->p[i].start; j < A->p[i].end; j++)
        {
            AT->x[count[A->i[j]]] = A->x[j];
            AT->i[count[A->i[j]]] = i;
            count[A->i[j]]++;
        }
    }

    return AT;
}

int calc_memory(int nnz, int n_rows, int extra_row_space, double memory_ratio)
{
    return (int) (nnz * memory_ratio) + n_rows * extra_row_space;
}

int calc_memory_row(int size, int extra_row_space, double memory_ratio)
{
    return (int) (size * memory_ratio) + extra_row_space;
}

void free_matrix(Matrix *A)
{
    if (A)
    {
        PS_FREE(A->i);
        PS_FREE(A->p);
        PS_FREE(A->x);
    }

    PS_FREE(A);
}

void remove_extra_space(Matrix *A, const int *row_sizes, const int *col_sizes,
                        bool remove_all, int *col_idxs_map)
{
    int i, j, start, end, len, row_alloc, curr, n_deleted_rows, col_count;
    int extra_row_space = (remove_all) ? 0 : EXTRA_ROW_SPACE;
    double extra_mem_ratio = (remove_all) ? 1.0 : EXTRA_MEMORY_RATIO;
    curr = 0;
    n_deleted_rows = 0;

    // --------------------------------------------------------------------------
    // loop through the rows and remove redundant space, including inactive
    // rows.
    // --------------------------------------------------------------------------
    for (i = 0; i < A->m; ++i)
    {
        if (row_sizes[i] == SIZE_INACTIVE_ROW)
        {
            n_deleted_rows++;
            continue;
        }

        start = A->p[i].start;
        end = A->p[i].end;
        len = end - start;
        row_alloc = calc_memory_row(len, extra_row_space, extra_mem_ratio);
        memmove(A->x + curr, A->x + start, ((size_t) (len)) * sizeof(double));
        memmove(A->i + curr, A->i + start, ((size_t) (len)) * sizeof(int));
        A->p[i - n_deleted_rows].start = curr;
        A->p[i - n_deleted_rows].end = curr + len;
        curr += row_alloc;
    }

    A->m -= n_deleted_rows;
    A->p[A->m].start = curr;
    A->p[A->m].end = curr;

    // shrink size
    A->x = (double *) ps_realloc(A->x, (size_t) MAX(curr, 1), sizeof(double));
    A->i = (int *) ps_realloc(A->i, (size_t) MAX(curr, 1), sizeof(int));
    A->p = (RowRange *) ps_realloc(A->p, (size_t) (A->m + 1), sizeof(RowRange));

    // -------------------------------------------------------------------------
    //                      compute new column indices
    // -------------------------------------------------------------------------
    col_count = 0;
    for (i = 0; i < A->n; ++i)
    {
        if (col_sizes[i] == SIZE_INACTIVE_COL)
        {
            col_idxs_map[i] = -1;
        }
        else
        {
            col_idxs_map[i] = (col_count++);
        }
    }
    A->n = col_count;

    // -------------------------------------------------------------------------
    //                        update column indices
    // -------------------------------------------------------------------------
    for (i = 0; i < A->m; ++i)
    {
        for (j = A->p[i].start; j < A->p[i].end; ++j)
        {
            A->i[j] = col_idxs_map[A->i[j]];
        }
    }
}

bool shift_row(Matrix *A, int row, int extra_space, int max_shift)
{
    int left, right, missing_space, remaining_shifts, left_shifts;
    int right_shifts, space_left, space_right, n_move_right, n_move_left;
    int next_start, next_end;
    bool shift_left;
    RowRange *row_r = A->p;
    left = row;
    right = row + 1;
    remaining_shifts = max_shift;
    left_shifts = 0;
    right_shifts = 0;
    missing_space = extra_space - (row_r[right].start - row_r[row].end);

    if (missing_space <= 0)
    {
        return true;
    }

    // ------------------------------------------------------------------------
    // compute the new start index for row 'row' and a lower bound on the start
    // index for the the first active row following row 'row'.
    // ------------------------------------------------------------------------
    while (missing_space > 0)
    {
        if (left == 0 && right == A->m)
        {
            return false;
        }

        // space_left is the number of steps we can shift 'left' row to the
        // left without overwriting row 'left - 1'.
        // space_right is the number of steps we can shift 'right' row to
        // the right without overwriting row 'right + 1'.
        assert(left >= 0 && right <= A->m);
        space_left = (left == 0) ? 0 : row_r[left].start - row_r[left - 1].end;
        space_right =
            (right == A->m) ? 0 : row_r[right + 1].start - row_r[right].end;
        assert(space_left >= 0 && space_right >= 0);

        // number of elements that must be moved left resp. right if we shift
        // in a certain direction
        n_move_right = row_r[right].end - row_r[right].start;
        n_move_left = row_r[left].end - row_r[left].start;

        // decide which direction to shift
        if (left == 0)
        {
            if (right != A->m && n_move_right <= remaining_shifts)
            {
                shift_left = false;
            }
            else
            {
                return false;
            }
        }
        else if (right == A->m)
        {
            if (left != 0 && n_move_left <= remaining_shifts)
            {
                shift_left = true;
            }
            else
            {
                return false;
            }
        }
        else if (n_move_left == 0)
        {
            shift_left = true;
        }
        else if (n_move_right == 0)
        {
            shift_left = false;
        }
        else if (n_move_left <= remaining_shifts &&
                 (space_left / (double) (n_move_left) >=
                  space_right / (double) (n_move_right)))
        {
            shift_left = true;
        }
        else if (n_move_right <= remaining_shifts)
        {
            shift_left = false;
        }
        else
        {
            return false;
        }

        assert(!(shift_left && left == 0) && !(!shift_left && right == A->m));

        if (shift_left)
        {
            left_shifts = MIN(missing_space, space_left);
            missing_space -= left_shifts;
            remaining_shifts -= n_move_left;
            left -= 1;
        }
        else
        {
            right_shifts = MIN(missing_space, space_right);
            missing_space -= right_shifts;
            remaining_shifts -= n_move_right;
            right += 1;
        }
    }
    assert(remaining_shifts >= 0);

    // ------------------------------------------------------------------------
    //                 execute total left shift
    // ------------------------------------------------------------------------
    next_start = row_r[left + 1].start - left_shifts;
    for (; left < row; left++)
    {
        int len = row_r[left + 1].end - row_r[left + 1].start;
        if (len > 0)
        {
            memmove(A->x + next_start, A->x + row_r[left + 1].start,
                    ((size_t) (len)) * sizeof(double));
            memmove(A->i + next_start, A->i + row_r[left + 1].start,
                    ((size_t) (len)) * sizeof(int));
        }
        row_r[left + 1].start = next_start;
        row_r[left + 1].end = next_start + len;
        next_start += len;
    }

    // ------------------------------------------------------------------------
    //                 execute total right shift
    // ------------------------------------------------------------------------
    next_end = row_r[right - 1].end + right_shifts;
    for (; right > row + 1; right--)
    {
        int len = row_r[right - 1].end - row_r[right - 1].start;
        if (len > 0)
        {
            memmove(A->x + next_end - len, A->x + row_r[right - 1].start,
                    ((size_t) (len)) * sizeof(double));
            memmove(A->i + next_end - len, A->i + row_r[right - 1].start,
                    ((size_t) (len)) * sizeof(int));
        }
        row_r[right - 1].start = next_end - len;
        row_r[right - 1].end = next_end;
        next_end = row_r[right - 1].start;
    }

    assert(row_r[row + 1].start - row_r[row].end == extra_space);
    return true;
}

void print_row_starts(const RowRange *row_ranges, size_t len)
{
    for (size_t i = 0; i < len; ++i)
    {
        printf("%d ", row_ranges[i].start);
    }
    printf("\n");
}

double insert_or_update_coeff(Matrix *A, int row, int col, double val, int *row_size)
{
    int i, start, end, insertion;
    double old_val = 0.0;
    start = A->p[row].start;
    end = A->p[row].end;
    insertion = end;

    // -----------------------------------------------------------------
    //             find where it should be inserted
    // -----------------------------------------------------------------
    for (i = start; i < end; ++i)
    {
        if (A->i[i] >= col)
        {
            insertion = i;
            break;
        }
    }

    // -----------------------------------------------------------------
    // Insert the new value if it is nonzero. If it exists or should be
    // inserted in the end, we don't need to shift values.
    // -----------------------------------------------------------------
    if (ABS(val) > ZERO_TOL)
    {
        // assert(!IS_ZERO_FEAS_TOL(val));
        if (insertion == end)
        {
            A->x[insertion] = val;
            A->i[insertion] = col;
            A->p[row].end += 1;
            A->nnz += 1;
            *row_size += 1;
        }
        else if (A->i[insertion] == col)
        {
            old_val = A->x[insertion];
            A->x[insertion] = val;
        }
        else
        {
            memmove(A->x + insertion + 1, A->x + insertion,
                    ((size_t) (end - insertion)) * sizeof(double));
            memmove(A->i + insertion + 1, A->i + insertion,
                    ((size_t) (end - insertion)) * sizeof(int));

            // insert new value
            A->x[insertion] = val;
            A->i[insertion] = col;
            A->p[row].end += 1;
            A->nnz += 1;
            *row_size += 1;
        }
    }
    // if the new value is zero, we just have to shift
    else
    {
        // we only expect that the new value is zero if the coefficient
        // already exists
        assert(A->i[insertion] == col);

        // we only have to shift values if the zero is not in the end
        if (insertion != end - 1)
        {
            memmove(A->x + insertion, A->x + insertion + 1,
                    ((size_t) (end - insertion - 1)) * sizeof(double));
            memmove(A->i + insertion, A->i + insertion + 1,
                    ((size_t) (end - insertion - 1)) * sizeof(int));
        }

        A->p[row].end -= 1;
        A->nnz -= 1;
        *row_size -= 1;
    }

    assert(A->p[row].end <= A->p[row + 1].start);
    return old_val;
}

void remove_coeff(RowView *row, int col)
{
    int shift = 0;
    int len = *row->len;
    for (int i = 0; i < len; ++i)
    {
        if (row->cols[i] == col)
        {
            shift = 1;
        }

        row->vals[i] = row->vals[i + shift];
        row->cols[i] = row->cols[i + shift];
    }

    assert(shift != 0);
    (*row->range).end -= 1;
    *row->len -= 1;
}

void count_rows(const Matrix *A, int *row_sizes)
{
    for (int i = 0; i < A->m; ++i)
    {
        row_sizes[i] = A->p[i].end - A->p[i].start;
    }
}

#ifdef TESTING
// Function to create a random CSR matrix
Matrix *random_matrix_new(int n_rows, int n_cols, double density)
{
    // allocate memory
    int n_alloc_nnz = (int) (density * n_rows * n_cols);
    double *Ax = (double *) ps_malloc((size_t) n_alloc_nnz, sizeof(double));
    int *Ai = (int *) ps_malloc((size_t) n_alloc_nnz, sizeof(int));
    int *Ap = (int *) ps_malloc((size_t) (n_rows + 1), sizeof(int));
    if (!Ax || !Ai || !Ap)
    {
        PS_FREE(Ax);
        PS_FREE(Ai);
        PS_FREE(Ap);
        return NULL;
    }

    // Initialize random number generator
    srand(1);

    int nnz_count = 0; // Counter for nonzero elements
    Ap[0] = 0;

    for (int i = 0; i < n_rows; ++i)
    {
        int row_nnz = 0;

        // Randomly determine the number of nonzeros in this row
        for (int j = 0; j < n_cols; ++j)
        {
            if ((double) rand() / RAND_MAX < density)
            {
                if (nnz_count >= n_alloc_nnz)
                {
                    break;
                }

                Ax[nnz_count] = ((double) (rand() - rand()) / RAND_MAX) * 20.0;
                Ai[nnz_count] = j;
                ++nnz_count;
                ++row_nnz;
            }
        }
        Ap[i + 1] = Ap[i] + row_nnz;
    }

    // create matrix in modified CSR format
    Matrix *A = matrix_new(Ax, Ai, Ap, n_rows, n_cols, nnz_count);
    PS_FREE(Ax);
    PS_FREE(Ai);
    PS_FREE(Ap);

    return A;
}

void replace_row_A(Matrix *A, int row, double ratio, double *new_vals, int *cols_new,
                   int new_len)
{
    int i, len, start, n_new_elements;
    len = A->p[row].end - A->p[row].start;
    n_new_elements = new_len - len;

    // potentially shift row to get extra space
    if (n_new_elements > 0)
    {
        assert(shift_row(A, row, n_new_elements, 2000));
    }

    // replace the row
    start = A->p[row].start;
    for (i = 0; i < new_len; ++i)
    {
        A->x[start + i] = ratio * new_vals[i];
        A->i[start + i] = cols_new[i];
    }
    A->p[row].end = A->p[row].start + new_len;
    assert(A->p[row].end <= A->p[row + 1].start);
}

#endif

1	/*
2	* Copyright 2025 Daniel Cederberg
3	*
4	* This file is part of the PSLP project (LP Presolver).
5	*
6	* Licensed under the Apache License, Version 2.0 (the "License");
7	* you may not use this file except in compliance with the License.
8	* You may obtain a copy of the License at
9	*
10	* http://www.apache.org/licenses/LICENSE-2.0
11	*
12	* Unless required by applicable law or agreed to in writing, software
13	* distributed under the License is distributed on an "AS IS" BASIS,
14	* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15	* See the License for the specific language governing permissions and
16	* limitations under the License.
17	*/
18
19	#include "Matrix.h"
20	#include "Debugger.h"
21	#include "Memory_wrapper.h"
22	#include "Numerics.h"
23	#include "RowColViews.h"
24	#include "glbopts.h"
25	#include "stdlib.h"
26	#include "string.h"
27
28	Matrix matrix_new(const double Ax, const int Ai, const int Ap, int n_rows,	29✔
29	int n_cols, int nnz)
30	{
31	DEBUG(ASSERT_NO_ZEROS_D(Ax, ((size_t) (nnz))));	29✔
32	Matrix *A = matrix_alloc(n_rows, n_cols, nnz);	29✔
33	RETURN_PTR_IF_NULL(A, NULL);	29✔
34	int offset, i, row_size, row_alloc;
35
36	offset = 0;	29✔
37	for (i = 0; i < n_rows; ++i)	168✔
38	{
39	A->p[i].start = Ap[i] + offset;	139✔
40	memcpy(A->x + A->p[i].start, Ax + Ap[i],	139✔
41	((size_t) (Ap[i + 1] - Ap[i])) * sizeof(double));	139✔
42
43	memcpy(A->i + A->p[i].start, Ai + Ap[i],	139✔
44	((size_t) (Ap[i + 1] - Ap[i])) * sizeof(int));	139✔
45
46	A->p[i].end = Ap[i + 1] + offset;	139✔
47	row_size = A->p[i].end - A->p[i].start;	139✔
48	row_alloc = calc_memory_row(row_size, EXTRA_ROW_SPACE, EXTRA_MEMORY_RATIO);	139✔
49	offset += row_alloc - row_size;	139✔
50	}
51
52	A->p[n_rows].start = Ap[n_rows] + offset;	29✔
53	A->p[n_rows].end = A->p[n_rows].start;	29✔
54
55	return A;	29✔
56	}
57
58	// needed for the transpose function
59	Matrix *matrix_alloc(int n_rows, int n_cols, int nnz)	371✔
60	{
61	Matrix A = (Matrix ) ps_malloc(1, sizeof(Matrix));	371✔
62	RETURN_PTR_IF_NULL(A, NULL);	371✔
63
64	A->m = n_rows;	371✔
65	A->n = n_cols;	371✔
66	A->nnz = nnz;	371✔
67	A->n_alloc = calc_memory(nnz, n_rows, EXTRA_ROW_SPACE, EXTRA_MEMORY_RATIO);	371✔
68
69	#ifdef TESTING
70	A->i = (int *) ps_calloc(((size_t) (A->n_alloc)), sizeof(int));	371✔
71	A->p = (RowRange *) ps_calloc(((size_t) (n_rows + 1)), sizeof(RowRange));	371✔
72	A->x = (double *) ps_calloc(((size_t) (A->n_alloc)), sizeof(double));	371✔
73	#else
74	A->i = (int *) ps_malloc((size_t) (A->n_alloc), sizeof(int));
75	A->p = (RowRange *) ps_malloc((size_t) (n_rows + 1), sizeof(RowRange));
76	A->x = (double *) ps_malloc((size_t) (A->n_alloc), sizeof(double));
77	#endif
78
79	if (!A->i \|\| !A->p \|\| !A->x)	371✔
80	{
81	free_matrix(A);	×
82	return NULL;	×
83	}
84
85	return A;	371✔
86	}
87
88	Matrix matrix_new_no_extra_space(const double Ax, const int Ai, const int Ap,	68✔
89	int n_rows, int n_cols, int nnz)
90	{
91	Matrix A = (Matrix ) ps_malloc(1, sizeof(Matrix));	68✔
92	RETURN_PTR_IF_NULL(A, NULL);	68✔
93
94	A->m = n_rows;	68✔
95	A->n = n_cols;	68✔
96	A->nnz = nnz;	68✔
97	A->n_alloc = nnz;	68✔
98	A->i = (int *) ps_malloc((size_t) A->n_alloc, sizeof(int));	68✔
99	A->p = (RowRange *) ps_malloc((size_t) (n_rows + 1), sizeof(RowRange));	68✔
100	A->x = (double *) ps_malloc((size_t) A->n_alloc, sizeof(double));	68✔
101
102	if (!A->i \|\| !A->p \|\| !A->x)	68✔
103	{
104	free_matrix(A);	×
105	return NULL;	×
106	}
107
108	memcpy(A->x, Ax, ((size_t) (nnz)) * sizeof(double));	68✔
109	memcpy(A->i, Ai, ((size_t) (nnz)) * sizeof(int));	68✔
110
111	for (int i = 0; i <= n_rows; ++i)	441✔
112	{
113	A->p[i].start = Ap[i];	373✔
114	A->p[i].end = Ap[i + 1];	373✔
115	}
116
117	return A;	68✔
118	}
119
120	Matrix transpose(const Matrix A, int *work_n_cols)	342✔
121	{
122	Matrix *AT = matrix_alloc(A->n, A->m, A->nnz);	342✔
123	RETURN_PTR_IF_NULL(AT, NULL);	342✔
124	int i, j, start;
125	int *count = work_n_cols;	342✔
126	memset(count, 0, ((size_t) (A->n)) * sizeof(int));	342✔
127
128	// -------------------------------------------------------------------
129	// compute nnz in each column of A
130	// -------------------------------------------------------------------
131	for (i = 0; i < A->m; ++i)	1,761✔
132	{
133	for (j = A->p[i].start; j < A->p[i].end; ++j)	6,018✔
134	{
135	count[A->i[j]]++;	4,599✔
136	}
137	}
138	// ------------------------------------------------------------------
139	// compute row pointers, taking the extra space into account
140	// ------------------------------------------------------------------
141	AT->p[0].start = 0;	342✔
142	for (i = 0; i < A->n; ++i)	1,993✔
143	{
144	start = AT->p[i].start;	1,651✔
145	AT->p[i].end = start + count[i];	1,651✔
146	AT->p[i + 1].start =	1,651✔
147	start + calc_memory_row(count[i], EXTRA_ROW_SPACE, EXTRA_MEMORY_RATIO);	1,651✔
148	count[i] = start;	1,651✔
149	}
150
151	AT->p[A->n].start = AT->n_alloc;	342✔
152	AT->p[A->n].end = AT->n_alloc;	342✔
153
154	// ------------------------------------------------------------------
155	// fill transposed matrix (this is a bottleneck)
156	// ------------------------------------------------------------------
157	for (i = 0; i < A->m; ++i)	1,761✔
158	{
159	for (j = A->p[i].start; j < A->p[i].end; j++)	6,018✔
160	{
161	AT->x[count[A->i[j]]] = A->x[j];	4,599✔
162	AT->i[count[A->i[j]]] = i;	4,599✔
163	count[A->i[j]]++;	4,599✔
164	}
165	}
166
167	return AT;	342✔
168	}
169
170	int calc_memory(int nnz, int n_rows, int extra_row_space, double memory_ratio)	371✔
171	{
172	return (int) (nnz * memory_ratio) + n_rows * extra_row_space;	371✔
173	}
174
175	int calc_memory_row(int size, int extra_row_space, double memory_ratio)	2,365✔
176	{
177	return (int) (size * memory_ratio) + extra_row_space;	2,365✔
178	}
179
180	void free_matrix(Matrix *A)	439✔
181	{
182	if (A)	439✔
183	{
184	PS_FREE(A->i);	439✔
185	PS_FREE(A->p);	439✔
186	PS_FREE(A->x);	439✔
187	}
188
189	PS_FREE(A);	439✔
190	}	439✔
191
192	void remove_extra_space(Matrix A, const int row_sizes, const int *col_sizes,	155✔
193	bool remove_all, int *col_idxs_map)
194	{
195	int i, j, start, end, len, row_alloc, curr, n_deleted_rows, col_count;
196	int extra_row_space = (remove_all) ? 0 : EXTRA_ROW_SPACE;	155✔
197	double extra_mem_ratio = (remove_all) ? 1.0 : EXTRA_MEMORY_RATIO;	155✔
198	curr = 0;	155✔
199	n_deleted_rows = 0;	155✔
200
201	// --------------------------------------------------------------------------
202	// loop through the rows and remove redundant space, including inactive
203	// rows.
204	// --------------------------------------------------------------------------
205	for (i = 0; i < A->m; ++i)	890✔
206	{
207	if (row_sizes[i] == SIZE_INACTIVE_ROW)	735✔
208	{
209	n_deleted_rows++;	160✔
210	continue;	160✔
211	}
212
213	start = A->p[i].start;	575✔
214	end = A->p[i].end;	575✔
215	len = end - start;	575✔
216	row_alloc = calc_memory_row(len, extra_row_space, extra_mem_ratio);	575✔
217	memmove(A->x + curr, A->x + start, ((size_t) (len)) * sizeof(double));	575✔
218	memmove(A->i + curr, A->i + start, ((size_t) (len)) * sizeof(int));	575✔
219	A->p[i - n_deleted_rows].start = curr;	575✔
220	A->p[i - n_deleted_rows].end = curr + len;	575✔
221	curr += row_alloc;	575✔
222	}
223
224	A->m -= n_deleted_rows;	155✔
225	A->p[A->m].start = curr;	155✔
226	A->p[A->m].end = curr;	155✔
227
228	// shrink size
229	A->x = (double *) ps_realloc(A->x, (size_t) MAX(curr, 1), sizeof(double));	155✔
230	A->i = (int *) ps_realloc(A->i, (size_t) MAX(curr, 1), sizeof(int));	155✔
231	A->p = (RowRange *) ps_realloc(A->p, (size_t) (A->m + 1), sizeof(RowRange));	155✔
232
233	// -------------------------------------------------------------------------
234	// compute new column indices
235	// -------------------------------------------------------------------------
236	col_count = 0;	155✔
237	for (i = 0; i < A->n; ++i)	909✔
238	{
239	if (col_sizes[i] == SIZE_INACTIVE_COL)	754✔
240	{
241	col_idxs_map[i] = -1;	167✔
242	}
243	else
244	{
245	col_idxs_map[i] = (col_count++);	587✔
246	}
247	}
248	A->n = col_count;	155✔
249
250	// -------------------------------------------------------------------------
251	// update column indices
252	// -------------------------------------------------------------------------
253	for (i = 0; i < A->m; ++i)	730✔
254	{
255	for (j = A->p[i].start; j < A->p[i].end; ++j)	2,256✔
256	{
257	A->i[j] = col_idxs_map[A->i[j]];	1,681✔
258	}
259	}
260	}	155✔
261
262	bool shift_row(Matrix *A, int row, int extra_space, int max_shift)	20✔
263	{
264	int left, right, missing_space, remaining_shifts, left_shifts;
265	int right_shifts, space_left, space_right, n_move_right, n_move_left;
266	int next_start, next_end;
267	bool shift_left;
268	RowRange *row_r = A->p;	20✔
269	left = row;	20✔
270	right = row + 1;	20✔
271	remaining_shifts = max_shift;	20✔
272	left_shifts = 0;	20✔
273	right_shifts = 0;	20✔
274	missing_space = extra_space - (row_r[right].start - row_r[row].end);	20✔
275
276	if (missing_space <= 0)	20✔
277	{
UNCOV 278	return true;	×
279	}
280
281	// ------------------------------------------------------------------------
282	// compute the new start index for row 'row' and a lower bound on the start
283	// index for the the first active row following row 'row'.
284	// ------------------------------------------------------------------------
285	while (missing_space > 0)	52✔
286	{
287	if (left == 0 && right == A->m)	39✔
288	{
289	return false;	×
290	}
291
292	// space_left is the number of steps we can shift 'left' row to the
293	// left without overwriting row 'left - 1'.
294	// space_right is the number of steps we can shift 'right' row to
295	// the right without overwriting row 'right + 1'.
296	assert(left >= 0 && right <= A->m);	39✔
297	space_left = (left == 0) ? 0 : row_r[left].start - row_r[left - 1].end;	39✔
298	space_right =	39✔
299	(right == A->m) ? 0 : row_r[right + 1].start - row_r[right].end;	39✔
300	assert(space_left >= 0 && space_right >= 0);	39✔
301
302	// number of elements that must be moved left resp. right if we shift
303	// in a certain direction
304	n_move_right = row_r[right].end - row_r[right].start;	39✔
305	n_move_left = row_r[left].end - row_r[left].start;	39✔
306
307	// decide which direction to shift
308	if (left == 0)	39✔
309	{
310	if (right != A->m && n_move_right <= remaining_shifts)	7✔
311	{
312	shift_left = false;	5✔
313	}
314	else
315	{
316	return false;	2✔
317	}
318	}
319	else if (right == A->m)	32✔
320	{
321	if (left != 0 && n_move_left <= remaining_shifts)	5✔
322	{
323	shift_left = true;	4✔
324	}
325	else
326	{
327	return false;	1✔
328	}
329	}
330	else if (n_move_left == 0)	27✔
331	{
332	shift_left = true;	×
333	}
334	else if (n_move_right == 0)	27✔
335	{
336	shift_left = false;	11✔
337	}
338	else if (n_move_left <= remaining_shifts &&	16✔
339	(space_left / (double) (n_move_left) >=	11✔
340	space_right / (double) (n_move_right)))	11✔
341	{
342	shift_left = true;	2✔
343	}
344	else if (n_move_right <= remaining_shifts)	14✔
345	{
346	shift_left = false;	10✔
347	}
348	else
349	{
350	return false;	4✔
351	}
352
353	assert(!(shift_left && left == 0) && !(!shift_left && right == A->m));	32✔
354
355	if (shift_left)	32✔
356	{
357	left_shifts = MIN(missing_space, space_left);	6✔
358	missing_space -= left_shifts;	6✔
359	remaining_shifts -= n_move_left;	6✔
360	left -= 1;	6✔
361	}
362	else
363	{
364	right_shifts = MIN(missing_space, space_right);	26✔
365	missing_space -= right_shifts;	26✔
366	remaining_shifts -= n_move_right;	26✔
367	right += 1;	26✔
368	}
369	}
370	assert(remaining_shifts >= 0);	13✔
371
372	// ------------------------------------------------------------------------
373	// execute total left shift
374	// ------------------------------------------------------------------------
375	next_start = row_r[left + 1].start - left_shifts;	13✔
376	for (; left < row; left++)	19✔
377	{
378	int len = row_r[left + 1].end - row_r[left + 1].start;	6✔
379	if (len > 0)	6✔
380	{
381	memmove(A->x + next_start, A->x + row_r[left + 1].start,	6✔
382	((size_t) (len)) * sizeof(double));	6✔
383	memmove(A->i + next_start, A->i + row_r[left + 1].start,	6✔
384	((size_t) (len)) * sizeof(int));	6✔
385	}
386	row_r[left + 1].start = next_start;	6✔
387	row_r[left + 1].end = next_start + len;	6✔
388	next_start += len;	6✔
389	}
390
391	// ------------------------------------------------------------------------
392	// execute total right shift
393	// ------------------------------------------------------------------------
394	next_end = row_r[right - 1].end + right_shifts;	13✔
395	for (; right > row + 1; right--)	38✔
396	{
397	int len = row_r[right - 1].end - row_r[right - 1].start;	25✔
398	if (len > 0)	25✔
399	{
400	memmove(A->x + next_end - len, A->x + row_r[right - 1].start,	13✔
401	((size_t) (len)) * sizeof(double));	13✔
402	memmove(A->i + next_end - len, A->i + row_r[right - 1].start,	13✔
403	((size_t) (len)) * sizeof(int));	13✔
404	}
405	row_r[right - 1].start = next_end - len;	25✔
406	row_r[right - 1].end = next_end;	25✔
407	next_end = row_r[right - 1].start;	25✔
408	}
409
410	assert(row_r[row + 1].start - row_r[row].end == extra_space);	13✔
411	return true;	13✔
412	}
413
414	void print_row_starts(const RowRange *row_ranges, size_t len)	×
415	{
416	for (size_t i = 0; i < len; ++i)	×
417	{
418	printf("%d ", row_ranges[i].start);	×
419	}
420	printf("\n");	×
421	}	×
422
423	double insert_or_update_coeff(Matrix A, int row, int col, double val, int row_size)	58✔
424	{
425	int i, start, end, insertion;
426	double old_val = 0.0;	58✔
427	start = A->p[row].start;	58✔
428	end = A->p[row].end;	58✔
429	insertion = end;	58✔
430
431	// -----------------------------------------------------------------
432	// find where it should be inserted
433	// -----------------------------------------------------------------
434	for (i = start; i < end; ++i)	126✔
435	{
436	if (A->i[i] >= col)	124✔
437	{
438	insertion = i;	56✔
439	break;	56✔
440	}
441	}
442
443	// -----------------------------------------------------------------
444	// Insert the new value if it is nonzero. If it exists or should be
445	// inserted in the end, we don't need to shift values.
446	// -----------------------------------------------------------------
447	if (ABS(val) > ZERO_TOL)	58✔
448	{
449	// assert(!IS_ZERO_FEAS_TOL(val));
450	if (insertion == end)	48✔
451	{
452	A->x[insertion] = val;	2✔
453	A->i[insertion] = col;	2✔
454	A->p[row].end += 1;	2✔
455	A->nnz += 1;	2✔
456	*row_size += 1;	2✔
457	}
458	else if (A->i[insertion] == col)	46✔
459	{
460	old_val = A->x[insertion];	18✔
461	A->x[insertion] = val;	18✔
462	}
463	else
464	{
465	memmove(A->x + insertion + 1, A->x + insertion,	28✔
466	((size_t) (end - insertion)) * sizeof(double));	28✔
467	memmove(A->i + insertion + 1, A->i + insertion,	28✔
468	((size_t) (end - insertion)) * sizeof(int));	28✔
469
470	// insert new value
471	A->x[insertion] = val;	28✔
472	A->i[insertion] = col;	28✔
473	A->p[row].end += 1;	28✔
474	A->nnz += 1;	28✔
475	*row_size += 1;	28✔
476	}
477	}
478	// if the new value is zero, we just have to shift
479	else
480	{
481	// we only expect that the new value is zero if the coefficient
482	// already exists
483	assert(A->i[insertion] == col);	10✔
484
485	// we only have to shift values if the zero is not in the end
486	if (insertion != end - 1)	10✔
487	{
488	memmove(A->x + insertion, A->x + insertion + 1,	6✔
489	((size_t) (end - insertion - 1)) * sizeof(double));	6✔
490	memmove(A->i + insertion, A->i + insertion + 1,	6✔
491	((size_t) (end - insertion - 1)) * sizeof(int));	6✔
492	}
493
494	A->p[row].end -= 1;	10✔
495	A->nnz -= 1;	10✔
496	*row_size -= 1;	10✔
497	}
498
499	assert(A->p[row].end <= A->p[row + 1].start);	58✔
500	return old_val;	58✔
501	}
502
503	void remove_coeff(RowView *row, int col)	31✔
504	{
505	int shift = 0;	31✔
506	int len = *row->len;	31✔
507	for (int i = 0; i < len; ++i)	175✔
508	{
509	if (row->cols[i] == col)	144✔
510	{
511	shift = 1;	31✔
512	}
513
514	row->vals[i] = row->vals[i + shift];	144✔
515	row->cols[i] = row->cols[i + shift];	144✔
516	}
517
518	assert(shift != 0);	31✔
519	(*row->range).end -= 1;	31✔
520	*row->len -= 1;	31✔
521	}	31✔
522
523	void count_rows(const Matrix A, int row_sizes)	136✔
524	{
525	for (int i = 0; i < A->m; ++i)	780✔
526	{
527	row_sizes[i] = A->p[i].end - A->p[i].start;	644✔
528	}
529	}	136✔
530
531	#ifdef TESTING
532	// Function to create a random CSR matrix
UNCOV 533	Matrix *random_matrix_new(int n_rows, int n_cols, double density)	×
534	{
535	// allocate memory
UNCOV 536	int n_alloc_nnz = (int) (density * n_rows * n_cols);	×
NEW 537	double Ax = (double ) ps_malloc((size_t) n_alloc_nnz, sizeof(double));	×
NEW 538	int Ai = (int ) ps_malloc((size_t) n_alloc_nnz, sizeof(int));	×
NEW 539	int Ap = (int ) ps_malloc((size_t) (n_rows + 1), sizeof(int));	×
UNCOV 540	if (!Ax \|\| !Ai \|\| !Ap)	×
541	{
542	PS_FREE(Ax);	×
543	PS_FREE(Ai);	×
544	PS_FREE(Ap);	×
545	return NULL;	×
546	}
547
548	// Initialize random number generator
UNCOV 549	srand(1);	×
550
UNCOV 551	int nnz_count = 0; // Counter for nonzero elements	×
UNCOV 552	Ap[0] = 0;	×
553
UNCOV 554	for (int i = 0; i < n_rows; ++i)	×
555	{
UNCOV 556	int row_nnz = 0;	×
557
558	// Randomly determine the number of nonzeros in this row
UNCOV 559	for (int j = 0; j < n_cols; ++j)	×
560	{
UNCOV 561	if ((double) rand() / RAND_MAX < density)	×
562	{
UNCOV 563	if (nnz_count >= n_alloc_nnz)	×
564	{
UNCOV 565	break;	×
566	}
567
UNCOV 568	Ax[nnz_count] = ((double) (rand() - rand()) / RAND_MAX) * 20.0;	×
UNCOV 569	Ai[nnz_count] = j;	×
UNCOV 570	++nnz_count;	×
UNCOV 571	++row_nnz;	×
572	}
573	}
UNCOV 574	Ap[i + 1] = Ap[i] + row_nnz;	×
575	}
576
577	// create matrix in modified CSR format
UNCOV 578	Matrix *A = matrix_new(Ax, Ai, Ap, n_rows, n_cols, nnz_count);	×
UNCOV 579	PS_FREE(Ax);	×
UNCOV 580	PS_FREE(Ai);	×
UNCOV 581	PS_FREE(Ap);	×
582
UNCOV 583	return A;	×
584	}
585
UNCOV 586	void replace_row_A(Matrix A, int row, double ratio, double new_vals, int *cols_new,	×
587	int new_len)
588	{
589	int i, len, start, n_new_elements;
UNCOV 590	len = A->p[row].end - A->p[row].start;	×
UNCOV 591	n_new_elements = new_len - len;	×
592
593	// potentially shift row to get extra space
UNCOV 594	if (n_new_elements > 0)	×
595	{
UNCOV 596	assert(shift_row(A, row, n_new_elements, 2000));	×
597	}
598
599	// replace the row
UNCOV 600	start = A->p[row].start;	×
UNCOV 601	for (i = 0; i < new_len; ++i)	×
602	{
UNCOV 603	A->x[start + i] = ratio * new_vals[i];	×
UNCOV 604	A->i[start + i] = cols_new[i];	×
605	}
UNCOV 606	A->p[row].end = A->p[row].start + new_len;	×
UNCOV 607	assert(A->p[row].end <= A->p[row + 1].start);	×
UNCOV 608	}	×
609
610	#endif

dance858 / PSLP / 21378561283

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous