22558105036

Committed 02 Mar 2026 01:47AM UTC coverage: 89.515% (-0.02%) from 89.536%

Build # 22558105036

Build Type

Pull #44

github

Committed by

web-flow

Commit Message

Merge 3cf270376 into 33b8bb671

Pull Request Pull Request #44: remove explicit zeros

Coverage Stats

1276 of 1422 branches covered (89.73%)

Branch coverage included in aggregate %.

10 of 13 new or added lines in 1 file covered. (76.92%)

11 existing lines in 1 file now uncovered.

3966 of 4434 relevant lines covered (89.45%)

6096.82 hits per line

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

94.78

/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 "Binary_search.h"
#include "Debugger.h"
#include "Memory_wrapper.h"
#include "Numerics.h"
#include "PSLP_warnings.h"
#include "RowColViews.h"
#include "glbopts.h"
#include "stdlib.h"
#include "string.h"

/* forward declaration */
static inline void remove_explicit_zeros(Matrix *A);

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

    offset = 0;
    for (i = 0; i < n_rows; ++i)
    {
        A->p[i].start = Ap[i] + offset;
        len = (size_t) (Ap[i + 1] - Ap[i]);
        memcpy(A->x + A->p[i].start, Ax + Ap[i], len * sizeof(double));
        memcpy(A->i + A->p[i].start, Ai + Ap[i], len * 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(size_t n_rows, size_t n_cols, size_t 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(A->n_alloc, sizeof(int));
    A->p = (RowRange *) ps_calloc(n_rows + 1, sizeof(RowRange));
    A->x = (double *) ps_calloc(A->n_alloc, sizeof(double));
#else
    A->i = (int *) ps_malloc(A->n_alloc, sizeof(int));
    A->p = (RowRange *) ps_malloc(n_rows + 1, sizeof(RowRange));
    A->x = (double *) ps_malloc(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,
                                  size_t n_rows, size_t n_cols, size_t 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(A->n_alloc, sizeof(int));
    A->p = (RowRange *) ps_malloc(n_rows + 1, sizeof(RowRange));
    A->x = (double *) ps_malloc(A->n_alloc, sizeof(double));

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

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

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

    /* the presolver assumes that only nonzero entries are stored */
    remove_explicit_zeros(A);

    return A;
}

static inline void remove_explicit_zeros(Matrix *A)
{
    int i, j, shift;
    for (i = 0; i < A->m; ++i)
    {
        shift = 0;
        for (j = A->p[i].start; j < A->p[i].end; ++j)
        {
            if (A->x[j] == 0.0)
            {
                shift++;
            }
            else if (shift > 0)
            {
                A->x[j - shift] = A->x[j];
                A->i[j - shift] = A->i[j];
            }
        }
        A->p[i].end -= shift;
        A->nnz -= (size_t) shift;
    }
}

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, 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 = (int) AT->n_alloc;
    AT->p[A->n].end = (int) 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;
}

size_t calc_memory(size_t nnz, size_t n_rows, size_t extra_row_space,
                   double memory_ratio)
{
    /* disable conversion compiler warning*/
    PSLP_DIAG_PUSH();
    PSLP_DIAG_IGNORE_CONVERSION();

    /* intentional truncation */
    size_t result = (size_t) (nnz * memory_ratio) + n_rows * extra_row_space;

    /* enable conversion compiler warnings */
    PSLP_DIAG_POP();
    return result;
}

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, bool remove_all,
                        const int *col_idxs_map, size_t new_n_cols)
{
    int j, start, end, len, row_alloc, curr;
    int extra_row_space = (remove_all) ? 0 : EXTRA_ROW_SPACE;
    double extra_mem_ratio = (remove_all) ? 1.0 : EXTRA_MEMORY_RATIO;
    curr = 0;
    size_t i, n_deleted_rows;

    // --------------------------------------------------------------------------
    // loop through the rows and remove redundant space, including inactive
    // rows.
    // --------------------------------------------------------------------------
    n_deleted_rows = 0;
    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;
        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;
        row_alloc = calc_memory_row(len, extra_row_space, extra_mem_ratio);
        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));

    // -------------------------------------------------------------------------
    //                        update column indices
    // -------------------------------------------------------------------------
    A->n = new_n_cols;
    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)
{
    double old_val = 0.0;
    int start = A->p[row].start;
    int end = A->p[row].end;

    // -----------------------------------------------------------------
    //             find where it should be inserted
    // -----------------------------------------------------------------
    int rel_ins = sorted_lower_bound(A->i + start, end - start, col);
    int insertion = start + rel_ins;

    // -----------------------------------------------------------------
    // 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)
    {
        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
        {
            size_t len = (size_t) (end - insertion);
            memmove(A->x + insertion + 1, A->x + insertion, len * sizeof(double));
            memmove(A->i + insertion + 1, A->i + insertion, len * 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)
        {
            size_t len = (size_t) (end - insertion - 1);
            memmove(A->x + insertion, A->x + insertion + 1, len * sizeof(double));
            memmove(A->i + insertion, A->i + insertion + 1, len * 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(size_t n_rows, size_t n_cols, double density)
{
    // allocate memory

    /* disable conversion compiler warning*/
    PSLP_DIAG_PUSH();
    PSLP_DIAG_IGNORE_CONVERSION();
    /* intentional truncation */
    size_t n_alloc_nnz = (size_t) (density * n_rows * n_cols);

    /* enable conversion compiler warnings */
    PSLP_DIAG_POP();
    double *Ax = (double *) ps_malloc(n_alloc_nnz, sizeof(double));
    int *Ai = (int *) ps_malloc(n_alloc_nnz, sizeof(int));
    int *Ap = (int *) ps_malloc(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);

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

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

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

dance858 / PSLP / 22558105036

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