/* linalg/cholesky_common.c
*
* Copyright (C) 2016 Patrick Alken
*
* 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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
/*
cholesky_swap_rowcol()
Swap rows and columns i and j of symmetric matrix A, updating only the
lower triangle
Inputs: A - symmetric matrix A, stored in lower triangle
i - index i
j - index j
Return: success/error
*/
static int
cholesky_swap_rowcol(gsl_matrix * A, const size_t i, const size_t j)
{
if (i != j)
{
const size_t N = A->size1;
double *Aii, *Ajj;
size_t ii, jj, k;
/* organize so that ii < jj */
if (i < j)
{
ii = i;
jj = j;
}
else
{
ii = j;
jj = i;
}
/* swap subrows A(i,1:i-1) with A(j,1:i-1) */
for (k = 0; k < ii; ++k)
{
double *Aik = gsl_matrix_ptr(A, ii, k);
double *Ajk = gsl_matrix_ptr(A, jj, k);
SWAP(*Aik, *Ajk);
}
/* swap subrow A(j,i+1:j-1) with subcolumn A(i+1:j-1,i) */
for (k = ii + 1; k < jj; ++k)
{
double *Ajk = gsl_matrix_ptr(A, jj, k);
double *Aki = gsl_matrix_ptr(A, k, ii);
SWAP(*Ajk, *Aki);
}
/* swap subcolumns A(j+1:N,i) with A(j+1:N,j) */
for (k = jj + 1; k < N; ++k)
{
double *Aki = gsl_matrix_ptr(A, k, ii);
double *Akj = gsl_matrix_ptr(A, k, jj);
SWAP(*Aki, *Akj);
}
/* now swap diagonal elements A(i,i) and A(j,j) */
Aii = gsl_matrix_ptr(A, ii, ii);
Ajj = gsl_matrix_ptr(A, jj, jj);
SWAP(*Aii, *Ajj);
}
return GSL_SUCCESS;
}