|
Packit |
67cb25 |
/* linalg/cholesky_common.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2016 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
cholesky_swap_rowcol()
|
|
Packit |
67cb25 |
Swap rows and columns i and j of symmetric matrix A, updating only the
|
|
Packit |
67cb25 |
lower triangle
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - symmetric matrix A, stored in lower triangle
|
|
Packit |
67cb25 |
i - index i
|
|
Packit |
67cb25 |
j - index j
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
cholesky_swap_rowcol(gsl_matrix * A, const size_t i, const size_t j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (i != j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
double *Aii, *Ajj;
|
|
Packit |
67cb25 |
size_t ii, jj, k;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* organize so that ii < jj */
|
|
Packit |
67cb25 |
if (i < j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ii = i;
|
|
Packit |
67cb25 |
jj = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ii = j;
|
|
Packit |
67cb25 |
jj = i;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* swap subrows A(i,1:i-1) with A(j,1:i-1) */
|
|
Packit |
67cb25 |
for (k = 0; k < ii; ++k)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double *Aik = gsl_matrix_ptr(A, ii, k);
|
|
Packit |
67cb25 |
double *Ajk = gsl_matrix_ptr(A, jj, k);
|
|
Packit |
67cb25 |
SWAP(*Aik, *Ajk);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* swap subrow A(j,i+1:j-1) with subcolumn A(i+1:j-1,i) */
|
|
Packit |
67cb25 |
for (k = ii + 1; k < jj; ++k)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double *Ajk = gsl_matrix_ptr(A, jj, k);
|
|
Packit |
67cb25 |
double *Aki = gsl_matrix_ptr(A, k, ii);
|
|
Packit |
67cb25 |
SWAP(*Ajk, *Aki);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* swap subcolumns A(j+1:N,i) with A(j+1:N,j) */
|
|
Packit |
67cb25 |
for (k = jj + 1; k < N; ++k)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double *Aki = gsl_matrix_ptr(A, k, ii);
|
|
Packit |
67cb25 |
double *Akj = gsl_matrix_ptr(A, k, jj);
|
|
Packit |
67cb25 |
SWAP(*Aki, *Akj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now swap diagonal elements A(i,i) and A(j,j) */
|
|
Packit |
67cb25 |
Aii = gsl_matrix_ptr(A, ii, ii);
|
|
Packit |
67cb25 |
Ajj = gsl_matrix_ptr(A, jj, jj);
|
|
Packit |
67cb25 |
SWAP(*Aii, *Ajj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|