/* $Id: cpl_matrix.c,v 1.94 2007/11/27 15:12:31 cizzo Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: cizzo $ * $Date: 2007/11/27 15:12:31 $ * $Revision: 1.94 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif #include #include #include #include #include #include #include #include #include "cpl_errorstate.h" #include "cpl_tools.h" /** * @defgroup cpl_matrix Matrices * * This module provides functions to create, destroy and use a @em cpl_matrix. * The elements of a @em cpl_matrix with M rows and N columns are counted * from 0,0 to M-1,N-1. The matrix element 0,0 is the one at the upper left * corner of a matrix. The CPL matrix functions work properly only in the * case the matrices elements do not contain garbage (such as @c NaN or * infinity). * * @par Synopsis: * @code * #include * @endcode */ /**@{*/ #ifndef inline #define inline /* inline */ #endif #define dtiny(x, t) ((x) < 0.0 ? (x) >= -(t) : (x) <= (t)) /* * The cpl_matrix type: */ struct _cpl_matrix_ { int nc; int nr; double *m; }; /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ /* * Private methods: */ /*----------------------------------------------------------------------------*/ /** @brief Swap two (different) rows @param self The matrix @param row1 One row @param row2 Another row @see cpl_matrix_swap_rows() @note No error checking done here Optimized due to repeated usage in cpl_matrix_decomp_lu(): 1) Removed redundant error checking - this allows the compiler to make 'swap' a register variable. 2) Use a single index for loop control and adressing the two rows - this redues the instruction count. This can have an effect when the swap is not memory bound. */ /*----------------------------------------------------------------------------*/ inline static void swap_rows(cpl_matrix * self, int row1, int row2) { double swap; int ncol = self->nc; double * pos1 = self->m + ncol * row1; double * pos2 = self->m + ncol * row2; while (ncol--) { swap = pos1[ncol]; pos1[ncol] = pos2[ncol]; pos2[ncol] = swap; } } /* * @brief * Just read a matrix row into the passed buffer. * * @param matrix Matrix where to read the row from. * @param pos Row number. * @param row Allocated buffer. * * @return Nothing. * * This private function reads the content of a matrix row into an allocated * buffer. No checks are performed. */ static void read_row(cpl_matrix *matrix, double *row, int pos) { int i; for (i = 0, pos *= matrix->nc; i < matrix->nc; i++, pos++) row[i] = matrix->m[pos]; } /* * @brief * Just write into a matrix row the values in buffer. * * @param matrix Matrix where to write the buffer to. * @param pos Row number. * @param row Allocated buffer. * * @return Nothing. * * This private function writes the content of a buffer into a matrix row. * No checks are performed. */ static void write_row(cpl_matrix *matrix, double *row, int pos) { int i; for (i = 0, pos *= matrix->nc; i < matrix->nc; i++, pos++) matrix->m[pos] = *row++; } /* * @brief * Just swap values in buffer with values in a matrix row. * * @param matrix Matrix to access. * @param pos Row to access. * @param row Allocated buffer. * * @return Nothing. * * This private function exchanges the values in buffer with the values in * a chosen matrix row. No checks are performed. */ static void write_read_row(cpl_matrix *matrix, double *row, int pos) { int i; double swap; for (i = 0, pos *= matrix->nc; i < matrix->nc; i++, pos++) { swap = matrix->m[pos]; matrix->m[pos] = row[i]; row[i] = swap; } } /* * @brief * Just read a matrix column into the passed buffer. * * @param matrix Matrix where to read the column from. * @param pos Column number. * @param column Allocated buffer. * * @return Nothing. * * This private function reads the content of a matrix column into an * allocated buffer. No checks are performed. */ static void read_column(cpl_matrix *matrix, double *column, int pos) { int i; for (i = 0; i < matrix->nr; i++, pos += matrix->nc) column[i] = matrix->m[pos]; } /* * @brief * Just write into a matrix column the values in buffer. * * @param matrix Matrix where to write the buffer to. * @param pos Column number. * @param column Allocated buffer. * * @return Nothing. * * This private function writes the content of a buffer into a matrix column. * No checks are performed. */ static void write_column(cpl_matrix *matrix, double *column, int pos) { int i; int size = matrix->nr * matrix->nc; for (i = pos; i < size; i += matrix->nc) matrix->m[i] = *column++; } /* * @brief * Just swap values in buffer with values in a matrix column. * * @param matrix Matrix to access. * @param pos Column to access. * @param column Allocated buffer. * * @return Nothing. * * This private function exchanges the values in buffer with the values in * a chosen matrix column. No checks are performed. */ static void write_read_column(cpl_matrix *matrix, double *column, int pos) { int i; double swap; for (i = 0; i < matrix->nr; i++, pos += matrix->nc) { swap = matrix->m[pos]; matrix->m[pos] = column[i]; column[i] = swap; } } /* * Public methods: */ /** * @brief * Print a matrix. * * @param matrix The matrix to print * @param stream The output stream * * @return Nothing. * * This function is intended just for debugging. It just prints the * elements of a matrix, ordered in rows and columns, to the specified * stream or FILE pointer. If the specified stream is NULL, it is set * to @em stdout. The function used for printing is the standard C * @c fprintf(). */ void cpl_matrix_dump(const cpl_matrix *matrix, FILE *stream) { int i, j, k; if (stream == NULL) stream = stdout; fprintf(stream, " "); if (matrix == NULL) { fprintf(stream, "NULL matrix\n\n"); return; } for (i = 0; i < matrix->nc; i++) fprintf(stream, " %d", i); fprintf(stream, "\n"); for (k = 0, i = 0; i < matrix->nr; i++) { fprintf(stream, " %d ", i); for (j = 0; j < matrix->nc; j++, k++) fprintf(stream, " %6g", matrix->m[k]); /* fprintf(stream, " %+6.2f", matrix->m[k]); */ fprintf(stream, "\n"); } fprintf(stream, "\n"); } /** * @brief * Create a zero matrix of given size. * * @param rows Number of matrix rows. * @param columns Number of matrix columns. * * @return Pointer to new matrix, or @c NULL in case of error. * * @error * * * * * *
CPL_ERROR_ILLEGAL_INPUT * rows or columns are not positive numbers. *
* @enderror * * This function allocates and initialises to zero a matrix of given size. * To destroy this matrix the function @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_new(int rows, int columns) { cpl_matrix *matrix; if (rows < 1 || columns < 1) { cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); return NULL; } matrix = (cpl_matrix *)cpl_malloc(sizeof(cpl_matrix)); matrix->m = (double *)cpl_calloc(rows * columns, sizeof(double)); matrix->nr = rows; matrix->nc = columns; return matrix; } /** * @brief * Create a new matrix from existing data. * * @param data Existing data buffer. * @param rows Number of matrix rows. * @param columns Number of matrix columns. * * @return Pointer to new matrix, or @c NULL in case of error. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input data is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * rows or columns are not positive numbers. *
* @enderror * * This function creates a new matrix that will encapsulate the given * data. At any error condition, a @c NULL pointer would be returned. * Note that the size of the input data array is not checked in any way, * and it is expected to match the specified matrix sizes. The input * array is supposed to contain in sequence all the new @em cpl_matrix * rows. For instance, in the case of a 3x4 matrix, the input array * should contain 12 elements * @code * 0 1 2 3 4 5 6 7 8 9 10 11 * @endcode * that would correspond to the matrix elements * @code * 0 1 2 3 * 4 5 6 7 * 8 9 10 11 * @endcode * The data buffer is not copied, so it should not be deallocated while * the matrix is still in use: the function @c cpl_matrix_delete() would * take care of deallocating it. To avoid problems with the memory managment, * the specified data buffer should be allocated using the functions of * the @c cpl_memory module, and also statically allocated data should be * strictly avoided. If this were not the case, this matrix should not be * destroyed using @c cpl_matrix_delete(), but @c cpl_matrix_unwrap() * should be used instead; moreover, functions implying memory handling * (as @c cpl_matrix_set_size(), or @c cpl_matrix_delete_row() ) should not * be used. */ cpl_matrix *cpl_matrix_wrap(int rows, int columns, double *data) { cpl_matrix *matrix; if (rows < 1 || columns < 1) { cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); return NULL; } if (data == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return NULL; } matrix = (cpl_matrix *)cpl_malloc(sizeof(cpl_matrix)); matrix->m = data; matrix->nr = rows; matrix->nc = columns; return matrix; } /** * @brief * Delete a matrix. * * @param matrix Pointer to a matrix to be deleted. * * @return Nothing. * * This function frees all the memory associated to a matrix. * If @em matrix is @c NULL, nothing is done. */ void cpl_matrix_delete(cpl_matrix *matrix) { if (matrix) { cpl_free(matrix->m); cpl_free(matrix); } } /** * @brief * Delete a matrix, but not its data buffer. * * @param matrix Pointer to a matrix to be deleted. * * @return Pointer to the internal data buffer. * * This function deallocates all the memory associated to a matrix, * with the exception of its data buffer. This type of destructor * should be used on matrices created with the @c cpl_matrix_wrap() * constructor, if the data buffer specified then was not allocated * using the functions of the @c cpl_memory module. In such a case, the * data buffer should be deallocated separately. See the documentation * of the function @c cpl_matrix_wrap(). If @em matrix is @c NULL, * nothing is done, and a @c NULL pointer is returned. */ void *cpl_matrix_unwrap(cpl_matrix *matrix) { void *data = NULL; if (matrix) { data = matrix->m; cpl_free(matrix); } return data; } /** * @brief * Get the number of rows of a matrix. * * @param matrix Pointer to the matrix to examine. * * @return Number of matrix rows, or zero in case of failure. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * Determine the number of rows in a matrix. */ inline int cpl_matrix_get_nrow(const cpl_matrix *matrix) { if (matrix == NULL) { /* inline precludes use of cpl_ensure() due to use of __func__ */ (void) cpl_error_set("cpl_matrix_get_nrow", CPL_ERROR_NULL_INPUT); return 0; } return matrix->nr; } /** * @brief * Get the number of columns of a matrix. * * @param matrix Pointer to the matrix to examine. * * @return Number of matrix columns, or zero in case of failure. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * Determine the number of columns in a matrix. */ inline int cpl_matrix_get_ncol(const cpl_matrix *matrix) { if (matrix == NULL) { /* inline precludes use of cpl_ensure() due to use of __func__ */ (void) cpl_error_set("cpl_matrix_get_ncol", CPL_ERROR_NULL_INPUT); return 0; } return matrix->nc; } /** * @brief * Get the pointer to a matrix data buffer, or @c NULL in case of error. * * @param matrix Input matrix. * * @return Pointer to the matrix data buffer. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * A @em cpl_matrix object includes an array of values of type @em double. * This function returns a pointer to this internal array, whose first * element corresponds to the @em cpl_matrix element 0,0. The internal * array contains in sequence all the @em cpl_matrix rows. For instance, * in the case of a 3x4 matrix, the array elements * @code * 0 1 2 3 4 5 6 7 8 9 10 11 * @endcode * would correspond to the following matrix elements: * @code * 0 1 2 3 * 4 5 6 7 * 8 9 10 11 * @endcode * * @note * Use at your own risk: direct manipulation of matrix data rules out * any check performed by the matrix object interface, and may introduce * inconsistencies between the information maintained internally, and * the actual matrix data and structure. */ inline double *cpl_matrix_get_data(cpl_matrix *matrix) { if (matrix == NULL) { /* inline precludes use of cpl_ensure() due to use of __func__ */ (void) cpl_error_set("cpl_matrix_get_data", CPL_ERROR_NULL_INPUT); return NULL; } return matrix->m; } /** * @brief * Get the pointer to a matrix data buffer, or @c NULL in case of error. * * @param matrix Input matrix. * * @return Pointer to the matrix data buffer. * * @see cpl_matrix_get_data */ inline const double *cpl_matrix_get_data_const(const cpl_matrix *matrix) { if (matrix == NULL) { /* inline precludes use of cpl_ensure() due to use of __func__ */ (void) cpl_error_set("cpl_matrix_get_data_const", CPL_ERROR_NULL_INPUT); return NULL; } return matrix->m; } /** * @brief * Get the value of a matrix element. * * @param matrix Pointer to a matrix. * @param row Matrix element row. * @param column Matrix element column. * * @return Value of the specified matrix element, or 0.0 in case of error. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The accessed element is beyond the matrix boundaries. *
* @enderror * * Get the value of a matrix element. The matrix rows and columns are * counted from 0,0. */ double cpl_matrix_get(const cpl_matrix *matrix, int row, int column) { if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return 0.0; } if (row < 0 || row >= matrix->nr || column < 0 || column >= matrix->nc) { cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); return 0.0; } return matrix->m[column + row * matrix->nc]; } /** * @brief * Write a value to a matrix element. * * @param matrix Input matrix. * @param row Matrix element row. * @param column Matrix element column. * @param value Value to write. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The accessed element is beyond the matrix boundaries. *
* @enderror * * Write a value to a matrix element. The matrix rows and columns are * counted from 0,0. */ cpl_error_code cpl_matrix_set(cpl_matrix *matrix, int row, int column, double value) { if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (row < 0 || row >= matrix->nr || column < 0 || column >= matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); matrix->m[column + row * matrix->nc] = value; return CPL_ERROR_NONE; } /** * @brief * Make a copy of a matrix. * * @param matrix Matrix to be duplicated. * * @return Pointer to the new matrix, or @c NULL in case of error. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * A copy of the input matrix is created. To destroy the duplicated matrix * the function @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_duplicate(const cpl_matrix *matrix) { cpl_matrix *new_matrix = NULL; int size; if (matrix) { new_matrix = (cpl_matrix *)cpl_malloc(sizeof(cpl_matrix)); new_matrix->nr = matrix->nr; new_matrix->nc = matrix->nc; size = new_matrix->nr * new_matrix->nc; new_matrix->m = (double *)cpl_malloc(size * sizeof(double)); memcpy(new_matrix->m, matrix->m, size * sizeof(double)); } else cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return new_matrix; } /** * @brief * Extract a submatrix from a matrix. * * @param matrix Pointer to the input matrix. * @param start_row Matrix row where to begin extraction. * @param start_column Matrix column where to begin extraction. * @param step_row Step between extracted rows. * @param step_column Step between extracted columns. * @param nrows Number of rows to extract. * @param ncolumns Number of columns to extract. * * @return Pointer to the new matrix, or @c NULL in case of error. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The start position is outside the matrix boundaries. *
CPL_ERROR_ILLEGAL_INPUT * While nrows and ncolumns are greater than 1, * the specified steps are not positive. *
* @enderror * * The new matrix will include the @em nrows x @em ncolumns values read * from the input matrix elements starting from position (@em start_row, * @em start_column), in a grid of steps @em step_row and @em step_column. * If the extraction parameters exceed the input matrix boundaries, just * the overlap is returned, and this matrix would have sizes smaller than * @em nrows x @em ncolumns. To destroy the new matrix the function * @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_extract(const cpl_matrix *matrix, int start_row, int start_column, int step_row, int step_column, int nrows, int ncolumns) { cpl_matrix *new_matrix = NULL; int i, r, c; int size; int row, col; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return NULL; } if (start_row < 0 || start_row >= matrix->nr || start_column < 0 || start_column >= matrix->nc || nrows < 1 || ncolumns < 1) { cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); return NULL; } if (nrows == 1) /* Step irrelevant... */ step_row = 1; if (ncolumns == 1) /* Step irrelevant... */ step_column = 1; if (step_row < 1 || step_column < 1) { cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); return NULL; } if ((start_row + (nrows - 1) * step_row) >= matrix->nr) nrows = (matrix->nr - start_row - 1) / step_row + 1; if ((start_column + (ncolumns - 1) * step_column) >= matrix->nc) ncolumns = (matrix->nc - start_column - 1) / step_column + 1; size = nrows * ncolumns; new_matrix = (cpl_matrix *)cpl_malloc(sizeof(cpl_matrix)); new_matrix->m = (double *)cpl_malloc(size * sizeof(double)); new_matrix->nr = nrows; new_matrix->nc = ncolumns; for (i = 0, r = start_row, row = nrows; row--; r += step_row) for (c = r * matrix->nc + start_column, col = ncolumns; col--; c += step_column, i++) new_matrix->m[i] = matrix->m[c]; return new_matrix; } /** * @brief * Extract a matrix row. * * @param matrix Pointer to the matrix containing the row. * @param row Sequence number of row to copy. * * @return Pointer to new matrix, or @c NULL in case of error. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The row is outside the matrix boundaries. *
* @enderror * * If a MxN matrix is given in input, the extracted row is a new 1xN matrix. * The row number is counted from 0. To destroy the new matrix the function * @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_extract_row(const cpl_matrix *matrix, int row) { cpl_matrix *vector; vector = cpl_matrix_extract(matrix, row, 0, 1, 1, 1, matrix->nc); if (vector == NULL) cpl_error_set_where("cpl_matrix_extract_row"); return vector; } /** * @brief * Copy a matrix column. * * @param matrix Pointer to matrix containing the column. * @param column Sequence number of column to copy. * * @return Pointer to new matrix, or @c NULL in case of error. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The column is outside the matrix boundaries. *
* @enderror * * If a MxN matrix is given in input, the extracted row is a new Mx1 matrix. * The column number is counted from 0. To destroy the new matrix the * function @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_extract_column(const cpl_matrix *matrix, int column) { cpl_matrix *vector; vector = cpl_matrix_extract(matrix, 0, column, 1, 1, matrix->nr, 1); if (vector == NULL) cpl_error_set_where("cpl_matrix_extract_column"); return vector; } /** * @brief * Extract a matrix diagonal. * * @param matrix Pointer to the matrix containing the diagonal. * @param diagonal Sequence number of the diagonal to copy. * * @return Pointer to the new matrix, or @c NULL in case of error. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The diagonal is outside the matrix boundaries * (see description below). *
* @enderror * * If a MxN matrix is given in input, the extracted diagonal is a Mx1 * matrix if N >= M, or a 1xN matrix if N < M. The diagonal number is * counted from 0, corresponding to the matrix diagonal starting at * element (0,0). A square matrix has just one diagonal; if M != N, * the number of diagonals in the matrix is |M - N| + 1. To specify * a @em diagonal sequence number outside this range would set an * error condition, and a @c NULL pointer would be returned. To destroy * the new matrix the function @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_extract_diagonal(const cpl_matrix *matrix, int diagonal) { cpl_matrix *output; int diagonal_count; int diagonal_length; int i, r, c; double *m; double *om; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return NULL; } diagonal_count = abs(matrix->nr - matrix->nc) + 1; if (diagonal < 0 || diagonal >= diagonal_count) { cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); return NULL; } diagonal_length = matrix->nr < matrix->nc ? matrix->nr : matrix->nc; output = (cpl_matrix *)cpl_malloc(sizeof(cpl_matrix)); om = output->m = (double *)cpl_malloc(diagonal_length * sizeof(double)); if (matrix->nc < matrix->nr) { r = diagonal; c = 0; output->nr = 1; output->nc = diagonal_length; } else { r = 0; c = diagonal; output->nr = diagonal_length; output->nc = 1; } m = matrix->m + c + r * matrix->nc; for (i = 0; i < diagonal_length; i++, m += matrix->nc + 1) *om++ = *m; return output; } /** * @brief * Write the same value to all matrix elements. * * @param matrix Pointer to the matrix to access * @param value Value to write * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * Write the same value to all matrix elements. */ cpl_error_code cpl_matrix_fill(cpl_matrix *matrix, double value) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; while (size--) *m++ = value; return CPL_ERROR_NONE; } /** * @brief * Write the same value to a matrix row. * * @param matrix Matrix to access * @param value Value to write * @param row Sequence number of row to overwrite. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified row is outside the matrix boundaries. *
* @enderror * * Write the same value to a matrix row. Rows are counted starting from 0. */ cpl_error_code cpl_matrix_fill_row(cpl_matrix *matrix, double value, int row) { double *m; int count; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (row < 0 || row >= matrix->nr) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); m = matrix->m + row * matrix->nc; count = matrix->nc; while (count--) *m++ = value; return CPL_ERROR_NONE; } /** * @brief * Write the same value to a matrix column. * * @param matrix Pointer to the matrix to access * @param value Value to write * @param column Sequence number of column to overwrite * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified column is outside the matrix boundaries. *
* @enderror * * Write the same value to a matrix column. Columns are counted starting * from 0. */ cpl_error_code cpl_matrix_fill_column(cpl_matrix *matrix, double value, int column) { int count; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (column < 0 || column >= matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); count = matrix->nr; m = matrix->m + column; while (count--) { *m = value; m += matrix->nc; } return CPL_ERROR_NONE; } /** * @brief * Write the same value to a matrix diagonal. * * @param matrix Matrix to access * @param value Value to write * @param diagonal Sequence number of diagonal to overwrite * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified diagonal is outside the matrix * boundaries (see description below). *
* @enderror * * Write the same value to a matrix diagonal. The diagonal number is * counted from 0, corresponding to the diagonal starting at element * (0,0). A square matrix has just one diagonal; if M != N, the number * of diagonals in the matrix is |M - N| + 1, counted along the * longer side of the matrix. */ cpl_error_code cpl_matrix_fill_diagonal(cpl_matrix *matrix, double value, int diagonal) { int diagonal_count; int diagonal_length; int i, r, c; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); diagonal_count = abs(matrix->nr - matrix->nc) + 1; if (diagonal < 0 || diagonal >= diagonal_count) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); diagonal_length = matrix->nr < matrix->nc ? matrix->nr : matrix->nc; if (matrix->nc < matrix->nr) { r = diagonal; c = 0; } else { r = 0; c = diagonal; } m = matrix->m + c + r * matrix->nc; for (i = 0; i < diagonal_length; i++, m += matrix->nc + 1) *m = value; return CPL_ERROR_NONE; } /** * @brief * Write the values of a matrix into another matrix. * * @param matrix Pointer to matrix to be modified. * @param submatrix Pointer to matrix to get the values from. * @param row Position of row 0 of @em submatrix in @em matrix. * @param col Position of column 0 of @em submatrix in @em matrix. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix or submatrix are NULL * pointers. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * No overlap exists between the two matrices. *
* @enderror * * The values of @em submatrix are written to @em matrix starting at the * indicated row and column. There are no restrictions on the sizes of * @em submatrix: just the parts of @em submatrix overlapping @em matrix * are copied. There are no restrictions on @em row and @em col either, * that can also be negative. If the two matrices do not overlap, nothing * is done, but an error condition is set. */ cpl_error_code cpl_matrix_copy(cpl_matrix *matrix, const cpl_matrix *submatrix, int row, int col) { int endrow; int endcol; int subrow; int subcol; int r, c, sr, sc; double *m; double *sm; if (matrix == NULL || submatrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); endrow = row + submatrix->nr; endcol = col + submatrix->nc; /* * Check whether matrices overlap. */ if (row >= matrix->nr || endrow < 1 || col >= matrix->nc || endcol < 1) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); /* * Define the overlap region on both matrices reference system: */ if (row < 0) { subrow = -row; row = 0; } else subrow = 0; if (col < 0) { subcol = -col; col = 0; } else subcol = 0; if (endrow > matrix->nr) endrow = matrix->nr; if (endcol > matrix->nc) endcol = matrix->nc; for (r = row, sr = subrow; r < endrow; r++, sr++) { m = matrix->m + r * matrix->nc + col; sm = submatrix->m + sr * submatrix->nc + subcol; for (c = col, sc = subcol; c < endcol; c++, sc++) *m++ = *sm++; } return CPL_ERROR_NONE; } /** * @brief * Write the same value into a submatrix of a matrix. * * @param matrix Pointer to matrix to be modified. * @param value Value to write. * @param row Start row of matrix submatrix. * @param col Start column of matrix submatrix. * @param nrow Number of rows of matrix submatrix. * @param ncol Number of columns of matrix submatrix. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified start position is outside the matrix * boundaries. *
CPL_ERROR_ILLEGAL_INPUT * nrow or ncol are not positive. *
* @enderror * * The specified value is written to @em matrix starting at the indicated * row and column; @em nrow and @em ncol can exceed the input matrix * boundaries, just the range overlapping the @em matrix is used in that * case. */ cpl_error_code cpl_matrix_fill_window(cpl_matrix *matrix, double value, int row, int col, int nrow, int ncol) { int endrow; int endcol; int r, c; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (row < 0 || row >= matrix->nr || col < 0 || col >= matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); if (nrow < 1 || ncol < 1) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); endrow = row + nrow; if (endrow > matrix->nr) endrow = matrix->nr; endcol = col + ncol; if (endcol > matrix->nc) endcol = matrix->nc; for (r = row; r < endrow; r++) { m = matrix->m + r * matrix->nc + col; for (c = col; c < endcol; c++) *m++ = value; } return CPL_ERROR_NONE; } /** * @brief * Shift matrix elements. * * @param matrix Pointer to matrix to be modified. * @param rshift Shift in the vertical direction. * @param cshift Shift in the horizontal direction. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The performed shift operation is cyclical (toroidal), i.e., matrix * elements shifted out of one side of the matrix get shifted in from * its opposite side. There are no restrictions on the values of the * shift. Positive shifts are always in the direction of increasing * row/column indexes. * * @note * The pointer to the matrix data buffer may change, therefore * pointers previously retrieved by calling @c cpl_matrix_get_data() * should be discarded. */ cpl_error_code cpl_matrix_shift(cpl_matrix *matrix, int rshift, int cshift) { cpl_matrix *shifted; double *m; double *s; int error; int size; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix->nr * matrix->nc == 1) return CPL_ERROR_NONE; rshift = rshift % matrix->nr; cshift = cshift % matrix->nc; if (rshift == 0 && cshift == 0) return CPL_ERROR_NONE; if (rshift < 0) rshift += matrix->nr; if (cshift < 0) cshift += matrix->nc; /* * This is not very efficient, but for the moment it may suffice... */ shifted = cpl_matrix_new(matrix->nr, matrix->nc); error = 0; if (cpl_matrix_copy(shifted, matrix, rshift, cshift)) error = 1; if (rshift) if (cpl_matrix_copy(shifted, matrix, rshift - shifted->nr, cshift)) error = 1; if (cshift) if (cpl_matrix_copy(shifted, matrix, rshift, cshift - shifted->nc)) error = 1; if (rshift && cshift) if (cpl_matrix_copy(shifted, matrix, rshift - shifted->nr, cshift - shifted->nc)) error = 1; if (error) { cpl_matrix_delete(shifted); cpl_error_set_where(cpl_func); error = cpl_error_get_code(); } else { m = matrix->m; s = shifted->m; size = matrix->nr * matrix->nc; while (size--) *m++ = *s++; cpl_matrix_delete(shifted); } return error; } /** * @brief * Rounding to zero very small numbers in matrix. * * @param matrix Pointer to matrix to be chopped. * @param tolerance Max tolerated rounding to zero. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * After specific manipulations of a matrix some of its elements * may theoretically be expected to be zero (for instance, as a result * of multiplying a matrix by its inverse). However, because of * numerical noise, such elements may turn out not to be exactly * zero. With this function any very small number in the matrix is * turned to exactly zero. If the @em tolerance is zero or negative, * a default threshold of @c DBL_EPSILON is used. */ cpl_error_code cpl_matrix_threshold_small(cpl_matrix *matrix, double tolerance) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; if (tolerance <= 0.0) tolerance = DBL_EPSILON; while (size--) { if (dtiny(*m, tolerance)) *m = 0.0; m++; } return CPL_ERROR_NONE; } /** * @brief * Check for zero matrix. * * @param matrix Pointer to matrix to be checked. * @param tolerance Max tolerated rounding to zero. * * @return 1 in case of zero matrix, 0 otherwise. If a @c NULL pointer * is passed, -1 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * After specific manipulations of a matrix some of its elements * may theoretically be expected to be zero. However, because of * numerical noise, such elements may turn out not to be exactly * zero. In this specific case, if any of the matrix element is * not exactly zero, the matrix would not be classified as a null * matrix. A threshold may be specified to consider zero any number * that is close enough to zero. If the specified @em tolerance is * negative, a default of @c DBL_EPSILON is used. A zero tolerance * may also be specified. */ int cpl_matrix_is_zero(const cpl_matrix *matrix, double tolerance) { int size; double *m; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return -1; } size = matrix->nr * matrix->nc; m = matrix->m; if (tolerance < 0.0) tolerance = DBL_EPSILON; while (size--) { if (!dtiny(*m, tolerance)) return 0; m++; } return 1; } /** * @brief * Check if a matrix is diagonal. * * @param matrix Pointer to matrix to be checked. * @param tolerance Max tolerated rounding to zero. * * @return 1 in case of diagonal matrix, 0 otherwise. If a @c NULL pointer * is passed, or the input matrix is not square, -1 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * A threshold may be specified to consider zero any number that * is close enough to zero. If the specified @em tolerance is * negative, a default of @c DBL_EPSILON is used. A zero tolerance * may also be specified. No error is set if the input matrix is * not square. */ int cpl_matrix_is_diagonal(const cpl_matrix *matrix, double tolerance) { int size; int dist; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return -1; } if (matrix->nr != matrix->nc) return 0; size = matrix->nr * matrix->nr; dist = matrix->nr + 1; if (tolerance < 0.0) tolerance = DBL_EPSILON; while (size--) if (size % dist) if (!dtiny(matrix->m[size], tolerance)) return 0; return 1; } /** * @brief * Check for identity matrix. * * @param matrix Pointer to matrix to be checked. * @param tolerance Max tolerated rounding to zero, or to one. * * @return 1 in case of identity matrix, 0 otherwise. If a @c NULL pointer * is passed, or the input matrix is not square, -1 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * A threshold may be specified to consider zero any number that is * close enough to zero, and 1 any number that is close enough to 1. * If the specified @em tolerance is negative, a default of @c DBL_EPSILON * is used. A zero tolerance may also be specified. No error is set if the * input matrix is not square. */ int cpl_matrix_is_identity(const cpl_matrix *matrix, double tolerance) { double tiny; int i; int skip; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return -1; } if (matrix->nr != matrix->nc) return 0; if (!cpl_matrix_is_diagonal(matrix, tolerance)) return 0; if (tolerance < 0.0) tolerance = DBL_EPSILON; i = 0; skip = matrix->nc + 1; for (i = 0; i < matrix->nr; i += skip) { tiny = matrix->m[i] - 1.0; if (!dtiny(tiny, tolerance)) return 0; } return 1; } /** * @brief * Swap two matrix rows. * * @param matrix Pointer to matrix to be modified. * @param row1 One matrix row. * @param row2 Another matrix row. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * Any of the specified rows is outside the matrix boundaries. *
* @enderror * * The values of two given matrix rows are exchanged. Rows are counted * starting from 0. If the same row number is given twice, nothing is * done and no error is set. */ cpl_error_code cpl_matrix_swap_rows(cpl_matrix *matrix, int row1, int row2) { cpl_ensure_code(matrix != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(row1 >= 0, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(row1 < matrix->nr, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(row2 >= 0, CPL_ERROR_ACCESS_OUT_OF_RANGE); cpl_ensure_code(row2 < matrix->nr, CPL_ERROR_ACCESS_OUT_OF_RANGE); if (row1 != row2) swap_rows(matrix, row1, row2); return CPL_ERROR_NONE; } /** * @brief * Swap two matrix columns. * * @param matrix Pointer to matrix to be modified. * @param column1 One matrix column. * @param column2 Another matrix column. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * Any of the specified columns is outside the matrix * boundaries. *
* @enderror * * The values of two given matrix columns are exchanged. Columns are * counted starting from 0. If the same column number is given twice, * nothing is done and no error is set. */ cpl_error_code cpl_matrix_swap_columns(cpl_matrix *matrix, int column1, int column2) { double swap; int nrow; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (column1 < 0 || column1 >= matrix->nc || column2 < 0 || column2 >= matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); if (column1 == column2) return CPL_ERROR_NONE; nrow = matrix->nr; while (nrow--) { swap = matrix->m[column1]; matrix->m[column1] = matrix->m[column2]; matrix->m[column2] = swap; column1 += matrix->nc; column2 += matrix->nc; } return CPL_ERROR_NONE; } /** * @brief * Swap a matrix column with a matrix row. * * @param matrix Pointer to matrix to be modified. * @param row Matrix row. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified row is outside the matrix boundaries. *
CPL_ERROR_ILLEGAL_INPUT * The input matrix is not square. *
* @enderror * * The values of the indicated row are exchanged with the column having * the same sequence number. Rows and columns are counted starting from 0. */ cpl_error_code cpl_matrix_swap_rowcolumn(cpl_matrix *matrix, int row) { double swap; int i; int posr, posc; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix->nr != matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); if (row < 0 || row >= matrix->nr) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); posr = row; posc = row * matrix->nc; for (i = 0; i < matrix->nr; i++) { swap = matrix->m[posr]; matrix->m[posr] = matrix->m[posc]; matrix->m[posc] = swap; posr += matrix->nc; posc++; } return CPL_ERROR_NONE; } /** * @brief * Reverse order of rows in matrix. * * @param matrix Pointer to matrix to be reversed. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The order of the rows in the matrix is reversed in place. */ cpl_error_code cpl_matrix_flip_rows(cpl_matrix *matrix) { int i, j; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); for (i = 0, j = matrix->nr - 1; i < j; i++, j--) swap_rows(matrix, i, j); return CPL_ERROR_NONE; } /** * @brief * Reverse order of columns in matrix. * * @param matrix Pointer to matrix to be reversed. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The order of the columns in the matrix is reversed in place. */ cpl_error_code cpl_matrix_flip_columns(cpl_matrix *matrix) { int i, j; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); for (i = 0, j = matrix->nc - 1; i < j; i++, j--) cpl_matrix_swap_columns(matrix, i, j); return CPL_ERROR_NONE; } /** * @brief * Create transposed matrix. * * @param matrix Pointer to matrix to be transposed. * * @return Pointer to transposed matrix. If a @c NULL pointer is passed, * a @c NULL pointer is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The transposed of the input matrix is created. To destroy the new matrix * the function @c cpl_matrix_delete() should be used. */ cpl_matrix *cpl_matrix_transpose_create(const cpl_matrix *matrix) { cpl_matrix *transposed = NULL; int i, j; double *m; double *tm; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return NULL; } tm = (double *)cpl_malloc(matrix->nc * matrix->nr * sizeof(double)); transposed = cpl_matrix_wrap(matrix->nc, matrix->nr, tm); m = matrix->m; for (i = 0; i < matrix->nr; i++) for (j = 0, tm = transposed->m + i; j < matrix->nc; j++, tm += matrix->nr) *tm = *m++; return transposed; } /** * @brief * Sort matrix by rows. * * @param matrix Pointer to matrix to be sorted. * @param mode Sorting mode: 0, by absolute value, otherwise by value. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The matrix elements of the leftmost column are used as reference for * the row sorting, if there are identical the values of the second * column are considered, etc. Rows with the greater values go on top. * If @em mode is equal to zero, the rows are sorted according to their * absolute values (zeroes at bottom). */ cpl_error_code cpl_matrix_sort_rows(cpl_matrix *matrix, int mode) { double *row; double value1, value2; int *done; int *sort_pattern; int *p; int i, j; int keep; int reach; int start; int count; int pos; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix->nr == 1) return CPL_ERROR_NONE; /* * Allocate the arrays for computing the sort pattern (i.e., the * list of arrival positions of array values after sorting). */ sort_pattern = cpl_malloc(matrix->nr * sizeof(int)); p = cpl_malloc(matrix->nr * sizeof(int)); for (i = 0; i < matrix->nr; i++) sort_pattern[i] = p[i] = i; /* * Find sorting pattern examining matrix columns. */ j = 0; i = 1; reach = 1; while (i < matrix->nr) { value1 = matrix->m[p[i] * matrix->nc + j]; value2 = matrix->m[p[i-1] * matrix->nc + j]; if (mode == 0) { value1 = fabs(value1); value2 = fabs(value2); } if (value1 < value2) { i++; if (i > reach) reach = i; } else if (value1 > value2) { keep = p[i-1]; p[i-1] = p[i]; p[i] = keep; if (i > 1) i--; else i = reach + 1; } else { if (matrix->nc == 1) { i++; continue; } j++; if (j == matrix->nc) { j = 0; i++; if (i > reach) reach = i; } continue; } if (j) { j = 0; continue; } } /* * To obtain the sort pattern this is the trick: */ reach = 1; i = 1; while (i < matrix->nr) { if (p[sort_pattern[i]] > p[sort_pattern[i-1]]) { i++; if (i > reach) reach = i; } else { keep = sort_pattern[i-1]; sort_pattern[i-1] = sort_pattern[i]; sort_pattern[i] = keep; if (i > 1) i--; else i = reach + 1; } } cpl_free(p); /* * Finally, sort the rows. Starting from row zero, see where it * should be written, save the overwritten row and write it to * the next position, etc. When a chain returns to the starting * point, the next undone matrix row is the new starting point, * and the process continues till all rows are done. This is * just a way to avoid doubling memory usage, paying with some * more CPU. */ done = (int *)cpl_calloc(matrix->nr, sizeof(int)); row = (double *)cpl_malloc(matrix->nc * sizeof(double)); pos = 0; start = 1; count = matrix->nr; while (count--) { done[pos] = 1; if (pos == sort_pattern[pos]) { if (count) for (pos = 0; pos < matrix->nr; pos++) if (!done[pos]) { start = 1; break; } continue; } if (start) { start = 0; read_row(matrix, row, pos); } pos = sort_pattern[pos]; if (done[pos]) { write_row(matrix, row, pos); if (count) { for (pos = 0; pos < matrix->nr; pos++) { if (!done[pos]) { start = 1; break; } } } } else write_read_row(matrix, row, pos); } cpl_free(sort_pattern); cpl_free(done); cpl_free(row); return CPL_ERROR_NONE; } /** * @brief * Sort matrix by columns. * * @param matrix Pointer to matrix to be sorted. * @param mode Sorting mode: 0, by absolute value, otherwise by value. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The matrix elements of the top row are used as reference for * the column sorting, if there are identical the values of the second * row are considered, etc. Columns with the largest values go on the * right. If @em mode is equal to zero, the columns are sorted according * to their absolute values (zeroes at left). */ cpl_error_code cpl_matrix_sort_columns(cpl_matrix *matrix, int mode) { double *column; double value1, value2; int *done; int *sort_pattern; int *p; int i, j; int keep; int reach; int start; int count; int pos; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix->nc == 1) return CPL_ERROR_NONE; /* * Allocate the arrays for computing the sort pattern (i.e., the * list of arrival positions of array values after sorting). */ sort_pattern = cpl_malloc(matrix->nc * sizeof(int)); p = cpl_malloc(matrix->nc * sizeof(int)); for (i = 0; i < matrix->nc; i++) sort_pattern[i] = p[i] = i; /* * Find sorting pattern examining matrix rows. */ j = 0; i = 1; reach = 1; while (i < matrix->nc) { value1 = matrix->m[j * matrix->nc + p[i]]; value2 = matrix->m[j * matrix->nc + p[i-1]]; if (mode == 0) { value1 = fabs(value1); value2 = fabs(value2); } if (value1 > value2) { i++; if (i > reach) reach = i; } else if (value1 < value2) { keep = p[i-1]; p[i-1] = p[i]; p[i] = keep; if (i > 1) i--; else i = reach + 1; } else { if (matrix->nr == 1) { i++; continue; } j++; if (j == matrix->nr) { j = 0; i++; if (i > reach) reach = i; } continue; } if (j) { j = 0; continue; } } /* * To obtain the sort pattern this is the trick: */ reach = 1; i = 1; while (i < matrix->nc) { if (p[sort_pattern[i]] > p[sort_pattern[i-1]]) { i++; if (i > reach) reach = i; } else { keep = sort_pattern[i-1]; sort_pattern[i-1] = sort_pattern[i]; sort_pattern[i] = keep; if (i > 1) i--; else i = reach + 1; } } cpl_free(p); /* * Finally, sort the columns. Starting from columns zero, see where * it should be written, save the overwritten column and write it * to the next position, etc. When a chain returns to the starting * point, the next undone matrix column is the new starting point, * and the process continues till all columns are done. This is * just a way to avoid doubling memory usage, paying with some * more CPU. */ done = (int *)cpl_calloc(matrix->nc, sizeof(int)); column = (double *)cpl_malloc(matrix->nr * sizeof(double)); pos = 0; start = 1; count = matrix->nc; while (count--) { done[pos] = 1; if (pos == sort_pattern[pos]) { if (count) for (pos = 0; pos < matrix->nc; pos++) if (!done[pos]) { start = 1; break; } continue; } if (start) { start = 0; read_column(matrix, column, pos); } pos = sort_pattern[pos]; if (done[pos]) { write_column(matrix, column, pos); if (count) { for (pos = 0; pos < matrix->nc; pos++) { if (!done[pos]) { start = 1; break; } } } } else write_read_column(matrix, column, pos); } cpl_free(sort_pattern); cpl_free(done); cpl_free(column); return CPL_ERROR_NONE; } /** * @brief * Delete rows from a matrix. * * @param matrix Pointer to matrix to be modified. * @param start First row to delete. * @param count Number of rows to delete. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified start is outside the matrix boundaries. *
CPL_ERROR_ILLEGAL_INPUT * count is not positive. *
CPL_ERROR_ILLEGAL_OUTPUT * Attempt to delete all the rows of matrix. *
* @enderror * * A portion of the matrix data is physically removed. The pointer * to matrix data may change, therefore pointers previously retrieved * by calling @c cpl_matrix_get_data() should be discarded. The specified * segment can extend beyond the end of the matrix, but the attempt to * remove all matrix rows is flagged as an error because zero length * matrices are illegal. Rows are counted starting from 0. */ cpl_error_code cpl_matrix_erase_rows(cpl_matrix *matrix, int start, int count) { double *m1; double *m2; int size; int i; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (count < 1) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); if (start < 0 || start >= matrix->nr) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); if (start + count > matrix->nr) count = matrix->nr - start; if (count == matrix->nr) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT); size = matrix->nr * matrix->nc; m1 = matrix->m + start * matrix->nc; m2 = m1 + count * matrix->nc; for (i = (start + count) * matrix->nc; i < size; i++) *m1++ = *m2++; size = (matrix->nr - count) * matrix->nc * sizeof(double); matrix->m = cpl_realloc(matrix->m, size); matrix->nr -= count; return CPL_ERROR_NONE; } /** * @brief * Delete columns from a matrix. * * @param matrix Pointer to matrix to be modified. * @param start First column to delete. * @param count Number of columns to delete. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ACCESS_OUT_OF_RANGE * The specified start is outside the matrix boundaries. *
CPL_ERROR_ILLEGAL_INPUT * count is not positive. *
CPL_ERROR_ILLEGAL_OUTPUT * Attempt to delete all the columns of matrix. *
* @enderror * * A portion of the matrix data is physically removed. The pointer * to matrix data may change, therefore pointers previously retrieved * by calling @c cpl_matrix_get_data() should be discarded. The specified * segment can extend beyond the end of the matrix, but the attempt to * remove all matrix columns is flagged as an error because zero length * matrices are illegal. Columns are counted starting from 0. */ cpl_error_code cpl_matrix_erase_columns(cpl_matrix *matrix, int start, int count) { double *m1; double *m2; int size; int new_nc; int i, j; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (count < 1) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); if (start < 0 || start >= matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ACCESS_OUT_OF_RANGE); if (start + count > matrix->nc) count = matrix->nc - start; if (count == matrix->nc) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT); new_nc = matrix->nc - count; for (j = 0; j < matrix->nr; j++) { m1 = matrix->m + j * new_nc; m2 = matrix->m + j * matrix->nc; for (i = 0; i < start; i++) *m1++ = *m2++; m2 += count; for (i += count; i < matrix->nc; i++) *m1++ = *m2++; } size = (matrix->nc - count) * matrix->nr * sizeof(double); matrix->m = cpl_realloc(matrix->m, size); matrix->nc -= count; return CPL_ERROR_NONE; } /** * @brief * Reframe a matrix. * * @param matrix Pointer to matrix to be modified. * @param top Extra rows on top. * @param bottom Extra rows on bottom. * @param left Extra columns on left. * @param right Extra columns on right. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_OUTPUT * Attempt to shrink matrix to zero size (or less). *
* @enderror * * The input matrix is reframed according to specifications. Extra rows * and column on the sides might also be negative, as long as they are * compatible with the matrix sizes: the input matrix would be reduced * in size accordingly, but an attempt to remove all matrix columns * and/or rows is flagged as an error because zero length matrices are * illegal. The old matrix elements contained in the new shape are left * unchanged, and new matrix elements added by the reshaping are initialised * to zero. No reshaping (i.e., all the extra rows set to zero) would * not be flagged as an error. * * @note * The pointer to the matrix data buffer may change, therefore * pointers previously retrieved by calling @c cpl_matrix_get_data() * should be discarded. */ cpl_error_code cpl_matrix_resize(cpl_matrix *matrix, int top, int bottom, int left, int right) { cpl_matrix *resized; int nr, nc; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (top == 0 && bottom == 0 && left == 0 && right == 0) return CPL_ERROR_NONE; nr = matrix->nr + top + bottom; nc = matrix->nc + left + right; if (nr < 1 || nc < 1) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT); resized = cpl_matrix_new(nr, nc); cpl_matrix_copy(resized, matrix, top, left); cpl_free(matrix->m); matrix->m = cpl_matrix_unwrap(resized); matrix->nr = nr; matrix->nc = nc; return CPL_ERROR_NONE; } /** * @brief * Resize a matrix. * * @param matrix Pointer to matrix to be resized. * @param rows New number of rows. * @param columns New number of columns. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_OUTPUT * Attempt to shrink matrix to zero size (or less). *
* @enderror * * The input matrix is resized according to specifications. The old * matrix elements contained in the resized matrix are left unchanged, * and new matrix elements added by an increase of the matrix number of * rows and/or columns are initialised to zero. * * @note * The pointer to the matrix data buffer may change, therefore * pointers previously retrieved by calling @c cpl_matrix_get_data() * should be discarded. */ cpl_error_code cpl_matrix_set_size(cpl_matrix *matrix, int rows, int columns) { if (cpl_matrix_resize(matrix, 0, rows - matrix->nr, 0, columns - matrix->nc)) { cpl_error_set_where("cpl_matrix_set_size"); return cpl_error_get_code(); } return CPL_ERROR_NONE; } /** * @brief * Append a matrix to another. * * @param matrix1 Pointer to first matrix. * @param matrix2 Pointer to second matrix. * @param mode Matrices connected horizontally (0) or vertically (1). * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * mode is neither 0 nor 1. *
CPL_ERROR_INCOMPATIBLE_INPUT * Matrices cannot be joined as indicated by mode. *
* @enderror * * If @em mode is set to 0, the matrices must have the same number * of rows, and are connected horizontally with the first matrix * on the left. If @em mode is set to 1, the matrices must have the * same number of columns, and are connected vertically with the * first matrix on top. The first matrix is expanded to include the * values from the second matrix, while the second matrix is left * untouched. * * @note * The pointer to the first matrix data buffer may change, therefore * pointers previously retrieved by calling @c cpl_matrix_get_data() * should be discarded. */ cpl_error_code cpl_matrix_append(cpl_matrix *matrix1, const cpl_matrix *matrix2, int mode) { int old_size; if (matrix1 == NULL || matrix2 == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (mode != 0 && mode != 1) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); if (mode == 0) { if (matrix1->nr != matrix2->nr) return cpl_error_set(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT); old_size = matrix1->nc; cpl_matrix_resize(matrix1, 0, 0, 0, matrix2->nc); cpl_matrix_copy(matrix1, matrix2, 0, old_size); } if (mode == 1) { if (matrix1->nc != matrix2->nc) return cpl_error_set(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT); old_size = matrix1->nr; cpl_matrix_resize(matrix1, 0, matrix2->nr, 0, 0); cpl_matrix_copy(matrix1, matrix2, old_size, 0); } return CPL_ERROR_NONE; } /** * @brief * Add two matrices. * * @param matrix1 Pointer to first matrix. * @param matrix2 Pointer to second matrix. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
CPL_ERROR_INCOMPATIBLE_INPUT * The specified matrices do not have the same size. *
* @enderror * * Add two matrices element by element. The two matrices must have * identical sizes. The result is written to the first matrix. */ cpl_error_code cpl_matrix_add(cpl_matrix *matrix1, const cpl_matrix *matrix2) { int size; double *m1; double *m2; if (matrix1 == NULL || matrix2 == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix1->nr != matrix2->nr || matrix1->nc != matrix2->nc) return cpl_error_set(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT); size = matrix1->nr * matrix1->nc; m1 = matrix1->m; m2 = matrix2->m; cpl_tools_add_flops( size ); while (size--) *m1++ += *m2++; return CPL_ERROR_NONE; } /** * @brief * Subtract a matrix from another. * * @param matrix1 Pointer to first matrix. * @param matrix2 Pointer to second matrix. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
CPL_ERROR_INCOMPATIBLE_INPUT * The specified matrices do not have the same size. *
* @enderror * * Subtract the second matrix from the first one element by element. * The two matrices must have identical sizes. The result is written * to the first matrix. */ cpl_error_code cpl_matrix_subtract(cpl_matrix *matrix1, const cpl_matrix *matrix2) { int size; double *m1; double *m2; if (matrix1 == NULL || matrix2 == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix1->nr != matrix2->nr || matrix1->nc != matrix2->nc) return cpl_error_set(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT); size = matrix1->nr * matrix1->nc; m1 = matrix1->m; m2 = matrix2->m; cpl_tools_add_flops( size ); while (size--) *m1++ -= *m2++; return CPL_ERROR_NONE; } /** * @brief * Multiply two matrices element by element. * * @param matrix1 Pointer to first matrix. * @param matrix2 Pointer to second matrix. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
CPL_ERROR_INCOMPATIBLE_INPUT * The specified matrices do not have the same size. *
* @enderror * * Multiply the two matrices element by element. The two matrices must * have identical sizes. The result is written to the first matrix. * * @note * To obtain the rows-by-columns product between two matrices, * @c cpl_matrix_product_create() should be used instead. */ cpl_error_code cpl_matrix_multiply(cpl_matrix *matrix1, const cpl_matrix *matrix2) { int size; double *m1; double *m2; if (matrix1 == NULL || matrix2 == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix1->nr != matrix2->nr || matrix1->nc != matrix2->nc) return cpl_error_set(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT); size = matrix1->nr * matrix1->nc; m1 = matrix1->m; m2 = matrix2->m; cpl_tools_add_flops( size ); while (size--) *m1++ *= *m2++; return CPL_ERROR_NONE; } /** * @brief * Divide a matrix by another element by element. * * @param matrix1 Pointer to first matrix. * @param matrix2 Pointer to second matrix. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
CPL_ERROR_INCOMPATIBLE_INPUT * The specified matrices do not have the same size. *
* @enderror * * Divide each element of the first matrix by the corresponding * element of the second one. The two matrices must have the same * number of rows and columns. The result is written to the first * matrix. No check is made against a division by zero. */ cpl_error_code cpl_matrix_divide(cpl_matrix *matrix1, const cpl_matrix *matrix2) { int size; double *m1; double *m2; if (matrix1 == NULL || matrix2 == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (matrix1->nr != matrix2->nr || matrix1->nc != matrix2->nc) return cpl_error_set(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT); size = matrix1->nr * matrix1->nc; m1 = matrix1->m; m2 = matrix2->m; cpl_tools_add_flops( size ); while (size--) *m1++ /= *m2++; return CPL_ERROR_NONE; } /** * @brief * Add a scalar to a matrix. * * @param matrix Pointer to matrix. * @param value Value to add. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * Add the same value to each matrix element. */ cpl_error_code cpl_matrix_add_scalar(cpl_matrix *matrix, double value) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( size ); while (size--) *m++ += value; return CPL_ERROR_NONE; } /** * @brief * Subtract a scalar to a matrix. * * @param matrix Pointer to matrix. * @param value Value to subtract. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * Subtract the same value to each matrix element. */ cpl_error_code cpl_matrix_subtract_scalar(cpl_matrix *matrix, double value) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( size ); while (size--) *m++ -= value; return CPL_ERROR_NONE; } /** * @brief * Multiply a matrix by a scalar. * * @param matrix Pointer to matrix. * @param value Multiplication factor. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * Multiply each matrix element by the same factor. */ cpl_error_code cpl_matrix_multiply_scalar(cpl_matrix *matrix, double value) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( size ); while (size--) *m++ *= value; return CPL_ERROR_NONE; } /** * @brief * Divide a matrix by a scalar. * * @param matrix Pointer to matrix. * @param value Divisor. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_DIVISION_BY_ZERO * The input value is 0.0. *
* @enderror * * Divide each matrix element by the same value. */ cpl_error_code cpl_matrix_divide_scalar(cpl_matrix *matrix, double value) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (value == 0.0) return cpl_error_set(cpl_func, CPL_ERROR_DIVISION_BY_ZERO); size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( size ); while (size--) *m++ /= value; return CPL_ERROR_NONE; } /** * @brief * Compute the logarithm of matrix elements. * * @param matrix Pointer to matrix. * @param base Logarithm base. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * The input base, or any matrix element, is * not positive. *
* @enderror * * Each matrix element is replaced by its logarithm in the specified base. * The base and all matrix elements must be positive. If this is not the * case, the matrix would not be modified. */ cpl_error_code cpl_matrix_logarithm(cpl_matrix *matrix, double base) { int size; double *m; double logbase; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (base <= 0.0) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; while (size--) { if (*m <= 0.0) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); m++; } logbase = log(base); size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( 1 + 2 * size ); while (size--) { *m = log(*m) / logbase; m++; } return CPL_ERROR_NONE; } /** * @brief * Compute the exponential of matrix elements. * * @param matrix Target matrix. * @param base Exponential base. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * The input base is not positive. *
* @enderror * * Each matrix element is replaced by its exponential in the specified base. * The base must be positive. */ cpl_error_code cpl_matrix_exponential(cpl_matrix *matrix, double base) { int size; double *m; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); if (base <= 0.0) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( size ); while (size--) { *m = pow(base, *m); m++; } return CPL_ERROR_NONE; } /** * @brief * Compute a power of matrix elements. * * @param matrix Pointer to matrix. * @param exponent Constant exponent. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * Any matrix element is not compatible with the * specified exponent (see description below). *
* @enderror * * Each matrix element is replaced by its power to the specified * exponent. If the specified exponent is not negative, all matrix * elements must be not negative; if the specified exponent is * negative, all matrix elements must be positive; otherwise, an * error condition is set and the matrix will be left unchanged. * If the exponent is exactly 0.5 the (faster) @c sqrt() will be * applied instead of @c pow(). If the exponent is zero, then any * (non negative) matrix element would be assigned the value 1.0. */ cpl_error_code cpl_matrix_power(cpl_matrix *matrix, double exponent) { int size; double *m; int negative = exponent < 0.0; if (matrix == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; m = matrix->m; if (negative) { while (size--) { if (*m <= 0.0) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); m++; } } else { while (size--) { if (*m < 0.0) return cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT); m++; } } if (negative) exponent = -exponent; size = matrix->nr * matrix->nc; m = matrix->m; cpl_tools_add_flops( size ); if (negative) { while (size--) { *m = 1 / pow(*m, exponent); m++; } } else { if (exponent == 0.5) { while (size--) { *m = sqrt(*m); m++; } } else { while (size--) { *m = pow(*m, exponent); m++; } } } return CPL_ERROR_NONE; } /** * @brief * Rows-by-columns product of two matrices. * * @param matrix1 Pointer to left side matrix. * @param matrix2 Pointer to right side matrix. * * @return Pointer to product matrix, or @c NULL in case of error. * * @error * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
CPL_ERROR_INCOMPATIBLE_INPUT * The number of columns of the first matrix is not equal to * the number of rows of the second matrix. *
* @enderror * * Rows-by-columns product of two matrices. The number of columns of the * first matrix must be equal to the number of rows of the second matrix. * To destroy the new matrix the function @c cpl_matrix_delete() should * be used. */ cpl_matrix *cpl_matrix_product_create(const cpl_matrix *matrix1, const cpl_matrix *matrix2) { cpl_matrix *self; cpl_ensure(matrix1 != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(matrix2 != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(matrix1->nc == matrix2->nr, CPL_ERROR_INCOMPATIBLE_INPUT, NULL); /* Create data-buffer without overhead of initializing */ self = cpl_matrix_wrap(matrix1->nr, matrix2->nc, cpl_malloc(matrix1->nr * matrix2->nc * sizeof(double))); (void) cpl_matrix_product(self, matrix1, matrix2); return self; } /* * @brief Replace a matrix by its LU-decomposition * @param self n X n non-singular matrix to decompose * @param perm n integers to be filled with the row permutations * @param psig On succes set to 1/-1 for an even/odd number of permutations * @return CPL_ERROR_NONE on success, or the relevant CPL error code * @note If the matrix is singular the elements of self become undefined * @see Golub & Van Loan, Matrix Computations, Algorithms 3.2.1 (Outer Product * Gaussian Elimination) and 3.4.1 (Gauss Elimination with Partial Pivoting). * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * An input pointer is NULL. *
CPL_ERROR_ILLEGAL_INPUT * self is not an n by n matrix. *
CPL_ERROR_SINGULAR_MATRIX * self is singular. *
* @enderror * */ cpl_error_code cpl_matrix_decomp_lu(cpl_matrix * self, int * perm, int * psig) { int n, i, j; const double * aread; double * awrite; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(perm != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(psig != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(self->nc == self->nr, CPL_ERROR_ILLEGAL_INPUT); n = self->nc; aread = awrite = self->m; /* No row permutations yet */ *psig = 1; for (i=0; i < n; i++) perm[i] = i; for (j = 0; j < n - 1; j++) { /* Find maximum in the j-th column */ double pivot = fabs(aread[j + j * n]); int ipivot = j; for (i = j + 1; i < n; i++) { const double aij = fabs(aread[j + i * n]); if (aij > pivot) { pivot = aij; ipivot = i; } } if (pivot <= 0.0) { return cpl_error_set_message_macro(cpl_func, CPL_ERROR_SINGULAR_MATRIX, __FILE__, __LINE__, "Pivot %d of %d is non-" "positive: %g", j, n, pivot); } if (ipivot > j) { /* The maximum is below the diagonal - swap rows */ const int iswap = perm[j]; perm[j] = perm[ipivot]; perm[ipivot] = iswap; *psig = - (*psig); swap_rows(self, j, ipivot); } pivot = aread[j + j * n]; for (i = j + 1; i < n; i++) { const double aij = aread[j + i * n] / pivot; int k; awrite[j + i * n] = aij; for (k = j + 1; k < n; k++) { const double ajk = aread[k + j * n]; const double aik = aread[k + i * n]; awrite[k + i * n] = aik - aij * ajk; } } } cpl_tools_add_flops( (2 * n * n * n) / 3); /* Check if A(n,n) is non-zero */ if (fabs(aread[n * (n-1) + (n-1)]) <= 0.0) { return cpl_error_set_message_macro(cpl_func, CPL_ERROR_SINGULAR_MATRIX, __FILE__, __LINE__, "The last of %d pivots is zero", n); } return CPL_ERROR_NONE; } /* * @brief Solve a LU-system * @param self n X n LU-matrix from cpl_matrix_decomp_lu() * @param rhs m right-hand-sides to be replaced by their solution * @param perm n integers filled with the row permutations, or NULL * @return CPL_ERROR_NONE on success, or the relevant CPL error code * @see cpl_matrix_decomp_lu() * * @error * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * An input pointer is NULL. *
CPL_ERROR_ILLEGAL_INPUT * self is not an n by n matrix. *
CPL_ERROR_INCOMPATIBLE_INPUT * The specified matrices do not have the same number of rows. *
CPL_ERROR_DIVISION_BY_ZERO * The main diagonal of U contains a zero. This error can only occur * if the LU-matrix does not come from a successful call to * cpl_matrix_decomp_lu(). *
* @enderror * */ cpl_error_code cpl_matrix_solve_lu(const cpl_matrix * self, cpl_matrix * rhs, const int * perm) { int n, i, j, k; double * x = NULL; /* Avoid false uninit warning */ const double * aread; const double * bread; double * bwrite; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(rhs != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(self->nc == self->nr, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(rhs->nr == self->nr, CPL_ERROR_INCOMPATIBLE_INPUT); n = self->nc; aread = self->m; bread = bwrite = rhs->m; if (perm != NULL) x = (double*) cpl_malloc(n * sizeof(double)); for (k=0; k < rhs->nc; k++) { if (perm != NULL) { /* Un-permute the rows in column k */ for (i = 0; i < n; i++) { x[i] = bread[rhs->nc * i + k]; } for (i = 0; i < n; i++) { bwrite[rhs->nc * i + k] = x[perm[i]]; } } /* Forward substitution in column k */ for (i = 1; i < n; i++) { double tmp = bread[rhs->nc * i + k]; for (j = 0; j < i; j++) { const double aij = aread[n * i + j]; tmp -= aij * bread[rhs->nc * j + k]; } bwrite[rhs->nc * i + k] = tmp; } /* Back substitution in column k */ for (i = n - 1; i >= 0; i--) { double tmp = bread[rhs->nc * i + k]; for (j = i + 1; j < n; j++) { const double aij = aread[n * i + j]; tmp -= aij * bread[rhs->nc * j + k]; } /* Check for a bug in the calling function */ if (aread[n * i + i] == 0.0) break; bwrite[rhs->nc * i + k] = tmp/aread[n * i + i]; } if (i >= 0) break; } if (perm != NULL) cpl_free(x); /* Flop count may be a bit too high if the below check fails */ cpl_tools_add_flops( k * 2 * n * n ); cpl_ensure_code(k == rhs->nc, CPL_ERROR_DIVISION_BY_ZERO); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Replace a matrix by its Cholesky-decomposition, L * transpose(L) = A @param self N by N symmetric positive-definite matrix to decompose @return CPL_ERROR_NONE on success, or the relevant CPL error code @note Only the upper triangle of self is read, L is written in lower triangle @note If the matrix is singular the elements of self become undefined @see Golub & Van Loan, Matrix Computations, Algorithm 4.2.1 (Cholesky: Gaxpy Version). @error
CPL_ERROR_NULL_INPUT An input pointer is NULL.
CPL_ERROR_ILLEGAL_INPUT self is not an n by n matrix.
CPL_ERROR_SINGULAR_MATRIX self is not symmetric, positive definite.
@enderror */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_matrix_decomp_chol(cpl_matrix * self) { int n, i, j; const double * aread; double * awrite; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); n = cpl_matrix_get_ncol(self); cpl_ensure_code(cpl_matrix_get_nrow(self) == n, CPL_ERROR_ILLEGAL_INPUT); aread = awrite = cpl_matrix_get_data(self); for (i = 0; i < n; i++) { for (j = i; j < n; j++) { double sub = aread[n * i + j]; int k; for (k = i-1; k >= 0; k--) { sub -= aread[n * i + k] * aread[n * j + k]; } if (j > i) { awrite[n * j + i] = sub/aread[n * i + i]; } else if (sub > 0.0) { awrite[n * i + i] = sqrt(sub); } else { return cpl_error_set_message_macro(cpl_func, CPL_ERROR_SINGULAR_MATRIX, __FILE__, __LINE__, "Pivot %d of %d is non-" "positive: %g", i, n, sub); } } } /* Dominant term: N^3 / 3 */ cpl_tools_add_flops(( n * ( 1 + n * ( 3 + 2 * n))) / 6); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Solve a L*transpose(L)-system @param self N by N L*transpose(L)-matrix from cpl_matrix_decomp_chol() @param rhs M right-hand-sides to be replaced by their solution @return CPL_ERROR_NONE on success, or the relevant CPL error code @see cpl_matrix_decomp_chol() @note Only the lower triangle of self is accessed @error
CPL_ERROR_NULL_INPUT An input pointer is NULL.
CPL_ERROR_ILLEGAL_INPUT self is not an n by n matrix.
CPL_ERROR_INCOMPATIBLE_INPUT The specified matrices do not have the same number of rows.
CPL_ERROR_DIVISION_BY_ZERO The main diagonal of L contains a zero. This error can only occur if the L*transpose(L)-matrix does not come from a successful call to cpl_matrix_decomp_chol().
@enderror */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_matrix_solve_chol(const cpl_matrix * self, cpl_matrix * rhs) { int n, i, j, k; int nrhs; const double * aread; const double * bread; double * bwrite; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(rhs != NULL, CPL_ERROR_NULL_INPUT); n = cpl_matrix_get_ncol(self); cpl_ensure_code(cpl_matrix_get_nrow(self) == n, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(cpl_matrix_get_nrow(rhs) == n, CPL_ERROR_INCOMPATIBLE_INPUT); nrhs = cpl_matrix_get_ncol(rhs); aread = cpl_matrix_get_data_const(self); bread = bwrite = cpl_matrix_get_data(rhs); for (k=0; k < nrhs; k++) { /* Forward substitution in column k */ for (i = 0; i < n; i++) { double sub = bread[nrhs * i + k]; for (j = i-1; j >= 0; j--) { sub -= aread[n * i + j] * bread[nrhs * j + k]; } cpl_ensure_code(aread[n * i + i] != 0.0, CPL_ERROR_DIVISION_BY_ZERO); bwrite[nrhs * i + k] = sub/aread[n * i + i]; } /* Back substitution in column k */ for (i = n-1; i >= 0; i--) { double sub = bread[nrhs * i + k]; for (j = i+1; j < n; j++) { sub -= aread[n * j + i] * bread[nrhs * j + k]; } bwrite[nrhs * i + k] = sub/aread[n * i + i]; } } cpl_tools_add_flops( 2 * n * n * nrhs); return CPL_ERROR_NONE; } /** * @brief * Compute the determinant of a matrix. * * @param matrix Pointer to n by n matrix. * * @return Matrix determinant. In case of error, 0.0 is returned. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * The input matrix is not square. *
CPL_ERROR_UNSPECIFIED * The input matrix is near-singular with a determinant so * close to zero that it cannot be represented by a double. *
* @enderror * * The input matrix must be a square matrix. In case of a 1x1 matrix, * the matrix single element value is returned. */ double cpl_matrix_get_determinant(const cpl_matrix *matrix) { cpl_matrix * self; int * perm; double determinant; int sig; int i; cpl_errorstate prevstate = cpl_errorstate_get(); cpl_error_code error; cpl_ensure(matrix != NULL, CPL_ERROR_NULL_INPUT, 0.0); cpl_ensure(matrix->nr == matrix->nc, CPL_ERROR_ILLEGAL_INPUT, 0.0); perm = (int*) cpl_malloc(matrix->nr * sizeof(int)); self = cpl_matrix_duplicate(matrix); error = cpl_matrix_decomp_lu(self, perm, &sig); cpl_free(perm); cpl_errorstate_set(prevstate); if (!error) { /* Use long double because the intermediate products may vary greatly */ cpl_long_double ldet = (cpl_long_double)sig; for (i = 0; i < matrix->nr; i++) { const double aii = cpl_matrix_get(self, i, i); assert(aii != 0.0); ldet *= (cpl_long_double)aii; } determinant = (double)ldet; #if 0 if (fabs(determinant) == 0.0) { /* could probably get -0.0 */ /* Try to use logarithm... */ determinant = 0.0; for (i = 0; i < matrix->nr; i++) { const double aii = cpl_matrix_get(self, i, i); if (aii < 0.0) sig = -sig; determinant += log(fabs(aii)); } cpl_msg_info(cpl_func,"log(fabs(det(A)))=%g (sign=%d, n=%d)", determinant, sig, matrix->nr); determinant = sig * exp(determinant); } #endif cpl_matrix_delete(self); cpl_ensure(determinant != 0.0, CPL_ERROR_UNSPECIFIED, 0.0); } else { cpl_matrix_delete(self); /* A failure here should not be possible */ cpl_ensure(error == CPL_ERROR_SINGULAR_MATRIX, error, 0.0); determinant = 0.0; } return determinant; } /** * @brief * Solution of a linear system. * * @param coeff A non-singular N by N matrix. * @param rhs A matrix containing one or more right-hand-sides. * @note rhs must have N rows and may contain more than one column, * which each represent an independent right-hand-side. * * @return A newly allocated solution matrix with the size as rhs, * or @c NULL on error. * * @error * * * * * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * coeff is not a square matrix. *
CPL_ERROR_INCOMPATIBLE_INPUT * coeff and rhs do not have the same number of rows. *
CPL_ERROR_SINGULAR_MATRIX * coeff is singular (to working precision). *
* @enderror * * Compute the solution of a system of N equations with N unknowns: * * coeff * X = rhs * * @em coeff must be an NxN matrix, and @em rhs a NxM matrix. * M greater than 1 means that multiple independent right-hand-sides * are solved for. * To destroy the solution matrix the function @c cpl_matrix_delete() * should be used. * */ cpl_matrix *cpl_matrix_solve(const cpl_matrix *coeff, const cpl_matrix *rhs) { cpl_matrix * lu; cpl_matrix * x; int * perm; int n, sig; cpl_error_code error; cpl_ensure(coeff != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(rhs != NULL, CPL_ERROR_NULL_INPUT, NULL); n = coeff->nc; cpl_ensure(coeff->nr == n, CPL_ERROR_ILLEGAL_INPUT, NULL); cpl_ensure(rhs->nr == n, CPL_ERROR_INCOMPATIBLE_INPUT, NULL); lu = cpl_matrix_duplicate(coeff); perm = (int*) cpl_malloc(n * sizeof(int)); if (cpl_matrix_decomp_lu(lu, perm, &sig)) { cpl_matrix_delete(lu); cpl_free(perm); cpl_ensure(0, cpl_error_get_code(), NULL); } x = cpl_matrix_duplicate(rhs); /* should not be able to fail at this point */ error = cpl_matrix_solve_lu(lu, x, perm); cpl_matrix_delete(lu); cpl_free(perm); if (error) { cpl_matrix_delete(x); cpl_ensure(0, cpl_error_get_code(), NULL); } return x; } /** * @brief * Find a matrix inverse. * * @param matrix Pointer to matrix to invert. * * @return Inverse matrix. In case of error a @c NULL is returned. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input is a NULL pointer. *
CPL_ERROR_ILLEGAL_INPUT * self is not an n by n matrix. *
CPL_ERROR_SINGULAR_MATRIX * matrix cannot be inverted. *
* @enderror * * The input must be a square matrix. To destroy the new matrix the * function @c cpl_matrix_delete() should be used. * * @note * When calling @c cpl_matrix_invert_create() with a nearly singular * matrix, it is possible to get a result containin NaN values without * any error code being set. */ cpl_matrix *cpl_matrix_invert_create(const cpl_matrix *matrix) { cpl_matrix * lu; cpl_matrix * inverse; int * perm; int n, i; cpl_error_code error; cpl_ensure(matrix != NULL, CPL_ERROR_NULL_INPUT, NULL); n = matrix->nc; cpl_ensure(matrix->nr == n, CPL_ERROR_ILLEGAL_INPUT, NULL); lu = cpl_matrix_duplicate(matrix); perm = (int*) cpl_malloc(n * sizeof(int)); if (cpl_matrix_decomp_lu(lu, perm, &i)) { cpl_matrix_delete(lu); cpl_free(perm); cpl_ensure(0, cpl_error_get_code(), NULL); } /* Should not be able to fail at this point */ inverse = cpl_matrix_new(n, n); /* Create an identity matrix with the rows permuted */ for (i = 0; i < n; i++) { (void)cpl_matrix_set(inverse, i, perm[i], 1.0); } cpl_free(perm); error = cpl_matrix_solve_lu(lu, inverse, NULL); cpl_matrix_delete(lu); if (error) { cpl_matrix_delete(inverse); cpl_ensure(0, cpl_error_get_code(), NULL); } return inverse; } /** * @brief * Solution of overdetermined linear equations in a least squares sense. * * @param coeff The N by M matrix of coefficients, where N >= M. * @param rhs An N by K matrix containing K right-hand-sides. * @note rhs may contain more than one column, * which each represent an independent right-hand-side. * * @return A newly allocated M by K solution matrix, * or @c NULL on error. * * @error * * * * * * * * * * * * * *
CPL_ERROR_NULL_INPUT * Any input is a NULL pointer. *
CPL_ERROR_INCOMPATIBLE_INPUT * coeff and rhs do not have the same number of rows. *
CPL_ERROR_SINGULAR_MATRIX * The matrix is (near) singular and a solution cannot be computed. *
* @enderror * * The following linear system of N equations and M unknowns * is given: * * @code * coeff * X = rhs * @endcode * * where @em coeff is the NxM matrix of the coefficients, @em X is the * MxK matrix of the unknowns, and @em rhs the NxK matrix containing * the K right hand side(s). * * The solution to the normal equations is known to be a least-squares * solution, i.e. the 2-norm of coeff * X - rhs is minimized by the * solution to * transpose(coeff) * coeff * X = transpose(coeff) * rhs. * * In the case that coeff is square (N is equal to M) it gives a faster * and more accurate result to use cpl_matrix_solve(). * * The solution matrix should be deallocated with the function * @c cpl_matrix_delete(). */ cpl_matrix *cpl_matrix_solve_normal(const cpl_matrix *coeff, const cpl_matrix *rhs) { cpl_matrix * solution; cpl_matrix * At; cpl_matrix * AtA; cpl_error_code error; cpl_ensure(coeff != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(rhs != NULL, CPL_ERROR_NULL_INPUT, NULL); cpl_ensure(rhs->nr == coeff->nr, CPL_ERROR_INCOMPATIBLE_INPUT, NULL); At = cpl_matrix_transpose_create(coeff); solution = cpl_matrix_product_create(At, rhs); AtA = cpl_matrix_product_normal_create(At); cpl_matrix_delete(At); error = cpl_matrix_solve_spd(AtA, solution); cpl_matrix_delete(AtA); if (error) { cpl_matrix_delete(solution); cpl_ensure(0, error, NULL); } return solution; } /** * @brief * Find the mean of all matrix elements. * * @param matrix Pointer to matrix. * * @return Mean. In case of error 0.0 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The mean value of all matrix elements is calculated. * * @note * This function works properly only in the case all the elements of * the input matrix do not contain garbage (such as @c NaN or infinity). */ double cpl_matrix_get_mean(const cpl_matrix *matrix) { int size; double mean = 0.0; double *pt; if (matrix == NULL) { cpl_error_set("cpl_matrix_get_mean", CPL_ERROR_NULL_INPUT); return 0.0; } size = matrix->nr * matrix->nc; if (size == 1) return matrix->m[0]; cpl_tools_add_flops( size ); pt = matrix->m; while (size--) mean += *pt++; return (mean / (matrix->nr * matrix->nc)); } /** * @brief * Find the standard deviation of matrix elements. * * @param matrix Pointer to matrix. * * @return Standard deviation. In case of error, or if a matrix is * 1x1, 0.0 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The standard deviation of all matrix elements is calculated. * * @note * This function works properly only in the case all the elements of * the input matrix do not contain garbage (such as @c NaN or infinity). */ double cpl_matrix_get_stdev(const cpl_matrix *matrix) { int size; double mean; double stdev = 0.0; double *pt; if (matrix == NULL) { cpl_error_set("cpl_matrix_get_stdev", CPL_ERROR_NULL_INPUT); return 0.0; } size = matrix->nr * matrix->nc; if (size == 1) return 0.0; mean = cpl_matrix_get_mean(matrix); cpl_tools_add_flops( 4 * size + 2); pt = matrix->m; while (size--) { stdev += (*pt - mean) * (*pt - mean); pt++; } return sqrt(stdev / (matrix->nr * matrix->nc - 1)); } /** * @brief * Find the median of matrix elements. * * @param matrix Pointer to matrix. * * @return Median. In case of error 0.0 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The median value of all matrix elements is calculated. */ double cpl_matrix_get_median(const cpl_matrix *matrix) { double *copybuf; int size; double median = 0.0; if (matrix == NULL) { cpl_error_set("cpl_matrix_get_median", CPL_ERROR_NULL_INPUT); return 0.0; } size = matrix->nr * matrix->nc; if (size == 1) return matrix->m[0]; copybuf = cpl_malloc(size * sizeof(double)); memcpy(copybuf, matrix->m, size * sizeof(double)); median = cpl_tools_get_median_double(copybuf, size); cpl_free(copybuf); return median; } /** * @brief * Find the minimum value of matrix elements. * * @param matrix Pointer to matrix. * * @return Minimum value. In case of error, 0.0 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The minimum value of all matrix elements is found. */ double cpl_matrix_get_min(const cpl_matrix *matrix) { int size; double min; double *m; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return 0.0; } size = matrix->nr * matrix->nc; if (size == 1) return matrix->m[0]; m = matrix->m; min = *m; while (--size) { m++; if (*m < min) min = *m; } return min; } /** * @brief * Find the maximum value of matrix elements. * * @param matrix Pointer to matrix. * * @return Maximum value. In case of error, 0.0 is returned. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The maximum value of all matrix elements is found. */ double cpl_matrix_get_max(const cpl_matrix *matrix) { int size; double max; double *m; if (matrix == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return 0.0; } size = matrix->nr * matrix->nc; if (size == 1) return matrix->m[0]; m = matrix->m; max = *m; while (--size) { m++; if (*m > max) max = *m; } return max; } /** * @brief * Find position of minimum value of matrix elements. * * @param matrix Input matrix. * @param row Returned row position of minimum. * @param column Returned column position of minimum. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The position of the minimum value of all matrix elements is found. * If more than one matrix element have a value corresponding to * the minimum, the lowest element row number is returned in @em row. * If more than one minimum matrix elements have the same row number, * the lowest element column number is returned in @em column. */ cpl_error_code cpl_matrix_get_minpos(const cpl_matrix *matrix, int *row, int *column) { int i, size; double min; double *m; if (matrix == NULL || row == NULL || column == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; *row = 0; *column = 0; if (size == 1) return CPL_ERROR_NONE; m = matrix->m; min = *m; for (i = 1; i < size; i++) { m++; if (*m < min) { min = *m; *row = i / matrix->nc; *column = i % matrix->nc; } } return CPL_ERROR_NONE; } /** * @brief * Find position of the maximum value of matrix elements. * * @param matrix Input matrix. * @param row Returned row position of maximum. * @param column Returned column position of maximum. * * @return @c CPL_ERROR_NONE on success. * * @error * * * * * *
CPL_ERROR_NULL_INPUT * The input matrix is a NULL pointer. *
* @enderror * * The position of the maximum value of all matrix elements is found. * If more than one matrix element have a value corresponding to * the maximum, the lowest element row number is returned in @em row. * If more than one maximum matrix elements have the same row number, * the lowest element column number is returned in @em column. */ cpl_error_code cpl_matrix_get_maxpos(const cpl_matrix *matrix, int *row, int *column) { int i, size; double max; double *m; if (matrix == NULL || row == NULL || column == NULL) return cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); size = matrix->nr * matrix->nc; *row = 0; *column = 0; if (size == 1) return CPL_ERROR_NONE; m = matrix->m; max = *m; for (i = 1; i < size; i++) { m++; if (*m > max) { max = *m; *row = i / matrix->nc; *column = i % matrix->nc; } } return CPL_ERROR_NONE; } /**@}*/ /*----------------------------------------------------------------------------*/ /** @brief Solution of Symmetric, Positive Definite system of linear equations. * @param self The N by N Symmetric, Positive Definite matrix. @param rhs An N by K matrix containing K right-hand-sides. @note rhs may contain more than one column, which each represent an independent right-hand-side. The solution is written in rhs. Only the upper triangular part of self is read, while the lower triangular part is modified. @return CPL_ERROR_NONE on success, or the relevant CPL error code * @error
CPL_ERROR_NULL_INPUT self or rhs is a NULL pointer.
CPL_ERROR_INCOMPATIBLE_INPUT self and rhs do not have the same number of rows.
CPL_ERROR_ILLEGAL_INPUT self is not an N by N matrix.
CPL_ERROR_SINGULAR_MATRIX self is (near) singular and a solution cannot be computed.
@enderror */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_matrix_solve_spd(cpl_matrix *self, cpl_matrix *rhs) { cpl_ensure_code(!cpl_matrix_decomp_chol(self), cpl_error_get_code()); cpl_ensure_code(!cpl_matrix_solve_chol(self, rhs), cpl_error_get_code()); return CPL_ERROR_NONE; } /** * @brief Create and compute A = B * transpose(B) * * @param self M x N Matrix * @return Pointer to created M x M product matrix, or @c NULL on error. * @note Only the upper triangle of A is computed, while the elements * below the main diagonal have undefined values. * @see cpl_matrix_product_create() * * @error * * * * * *
CPL_ERROR_NULL_INPUT * Any input matrix is a NULL pointer. *
* @enderror * * To destroy the new matrix the function @c cpl_matrix_delete() should * be used. */ cpl_matrix * cpl_matrix_product_normal_create(const cpl_matrix * self) { double sum; cpl_matrix * product; const double * ai = cpl_matrix_get_data_const(self); const double * aj; double * bwrite; const int m = cpl_matrix_get_nrow(self); const int n = cpl_matrix_get_ncol(self); int i, j, k; cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL); #if 0 /* Initialize all values to zero. This is done to avoid access of uninitilized memory, in case someone passes the matrix to for example cpl_matrix_dump(). */ product = cpl_matrix_new(m, m); bwrite = cpl_matrix_get_data(product); #else bwrite = (double *) cpl_malloc(m * m * sizeof(double)); product = cpl_matrix_wrap(m, m, bwrite); #endif /* The result at (i,j) is the dot-product of i'th and j'th row */ for (i = 0; i < m; i++, bwrite += m, ai += n) { aj = ai; /* aj points to first entry in j'th row */ for (j = i; j < m; j++, aj += n) { sum = 0.0; for (k = 0; k < n; k++) { sum += ai[k] * aj[k]; } bwrite[j] = sum; } } cpl_tools_add_flops(n * m * (m + 1)); return product; } /*----------------------------------------------------------------------------*/ /** @brief Fill a matrix with the product of A * B @param self The matrix to fill, is or else will be set to size M x N @param ma The matrix A, of size M x K @param mb The matrix B, of size K * N @return CPL_ERROR_NONE or the relevant CPL error code on error @note If upon entry, self has a size m *n, which differs from its size on return, then the result of preceeding calls to cpl_matrix_get_data() are invalid. @see cpl_matrix_product_create() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_matrix_product(cpl_matrix * self, const cpl_matrix * ma, const cpl_matrix * mb) { double sum; double * ds = cpl_matrix_get_data(self); const double * d1 = cpl_matrix_get_data_const(ma); const double * d2 = cpl_matrix_get_data_const(mb); const int nr = cpl_matrix_get_nrow(ma); const int nc = cpl_matrix_get_ncol(mb); const int nk = cpl_matrix_get_nrow(mb); int i, j, k; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(ma != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(mb != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(ma->nc == nk, CPL_ERROR_INCOMPATIBLE_INPUT); if (self->nr != nr || self->nc != nc) { if (self->nr * self->nc != nr * nc) { cpl_free(self->m); self->m = cpl_malloc(nr * nc * sizeof(double)); } self->nr = nr; self->nc = nc; } for (i = 0; i < nr; i++, d1 += nk) { for (j = 0; j < nc; j++) { sum = 0.0; for (k = 0; k < nk; k++) { sum += d1[k] * d2[nc * k + j]; } ds[nc * i + j] = sum; } } cpl_tools_add_flops( 2 * ma->nr * mb->nr * mb->nc ); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Fill a matrix with the product of A * B' @param self The matrix to fill, is or else will be set to size M x N @param ma The matrix A, of size M x K @param mb The matrix B, of size N x K @return CPL_ERROR_NONE or the relevant CPL error code on error @note The use of the transpose of B causes a more efficient memory access @note Changing the order of A and B is allowed, it transposes the result @see cpl_matrix_product_create() */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_matrix_product_transpose(cpl_matrix * self, const cpl_matrix * ma, const cpl_matrix * mb) { double sum; double * ds = cpl_matrix_get_data(self); const double * d1 = cpl_matrix_get_data_const(ma); const double * d2 = cpl_matrix_get_data_const(mb); const double * di; const int nr = cpl_matrix_get_nrow(ma); const int nc = cpl_matrix_get_nrow(mb); const int nk = cpl_matrix_get_ncol(mb); int i, j, k; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(ma != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(mb != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(ma->nc == nk, CPL_ERROR_INCOMPATIBLE_INPUT); if (self->nr != nr || self->nc != nc) { if (self->nr * self->nc != nr * nc) { cpl_free(self->m); self->m = cpl_malloc(nr * nc * sizeof(double)); } self->nr = nr; self->nc = nc; } for (i = 0; i < nr; i++, d1 += nk) { /* Since ma and mb are addressed in the same manner, they can use the same index, k */ di = d2; /* di points to first entry in i'th row */ for (j = 0; j < nc; j++, di += nk) { sum = 0.0; for (k = 0; k < nk; k++) { sum += d1[k] * di[k]; } ds[nc * i + j] = sum; } } cpl_tools_add_flops( 2 * ma->nr * mb->nr * mb->nc ); return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Solve a L*transpose(L)-system with a transposed Right Hand Side @param self N by N L*transpose(L)-matrix from cpl_matrix_decomp_chol() @param rhs M right-hand-sides to be replaced by their solution @return CPL_ERROR_NONE on success, or the relevant CPL error code @see cpl_matrix_solve_chol() @note Only the lower triangle of self is accessed @note The use of the transpose of rhs causes a more efficient memory access @error
CPL_ERROR_NULL_INPUT An input pointer is NULL.
CPL_ERROR_ILLEGAL_INPUT self is not an n by n matrix.
CPL_ERROR_INCOMPATIBLE_INPUT Selfs number of rows differs from rhs' number of columns.
CPL_ERROR_DIVISION_BY_ZERO The main diagonal of L contains a zero. This error can only occur if the L*transpose(L)-matrix does not come from a successful call to cpl_matrix_decomp_chol().
@enderror */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_matrix_solve_chol_transpose(const cpl_matrix * self, cpl_matrix * rhs) { const int n = cpl_matrix_get_ncol(self); const int nrhs = cpl_matrix_get_nrow(rhs); int i, j, k; const double * aread; const double * ai; double * bk; double sub; cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(rhs != NULL, CPL_ERROR_NULL_INPUT); cpl_ensure_code(self->nr == n, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(rhs->nc == n, CPL_ERROR_INCOMPATIBLE_INPUT); aread = self->m; /* bk points to first entry in k'th right hand side */ bk = rhs->m; for (k=0; k < nrhs; k++, bk += n) { /* Forward substitution in column k */ /* Since self and rhs are addressed in the same manner, they can use the same index, j */ ai = aread; /* ai points to first entry in i'th row */ for (i = 0; i < n; i++, ai += n) { sub = 0.0; for (j = 0; j < i; j++) { sub += ai[j] * bk[j]; } cpl_ensure_code(k > 0 || ai[j] != 0.0, CPL_ERROR_DIVISION_BY_ZERO); bk[j] = (bk[j] - sub) / ai[j]; } /* Back substitution in column k */ for (i = n-1; i >= 0; i--) { sub = bk[i]; for (j = i+1; j < n; j++) { sub -= aread[n * j + i] * bk[j]; } bk[i] = sub/aread[n * i + i]; } } cpl_tools_add_flops( 2 * n * n * nrhs); return CPL_ERROR_NONE; }