I have been trying to use Givens rotations via a combination of drotg and dlasr but have not been able to get it right. Can anyone please help me out. Please find below the self contained code of what I am trying to acomplish.
Many thanks in advance,
Best regards,
Giovanni
- Code: Select all
/*
* test.cc
*
* Created on: May 31, 2012
* Author: Giovanni Azua
*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl_lapack.h>
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
static const int ROWS = 3;
static const int COLS = 3;
static double matrix[ROWS][COLS];
static void print_matrix() {
int rows = ROWS;
int cols = COLS;
printf("matrix of dimensions %dx%d: \n", rows, cols);
for (int i=0; i < rows; ++i) {
for (int j=0; j < cols; ++j) {
printf("%f ", matrix[j][i]);
}
printf("\n");
}
}
int main(int argc, char** argv) {
// create a simple matrix in "column major" with 1's in the
// upper triangule and 2s in the band below the diagonal
// now I want to get rid of the 2s using Givens rotations
matrix[0][0] = 1; matrix[1][0] = 1; matrix[2][0] = 1;
matrix[0][1] = 2; matrix[1][1] = 1; matrix[2][1] = 1;
matrix[0][2] = 0; matrix[1][2] = 2; matrix[2][2] = 1;
print_matrix();
double c[COLS];
double s[COLS];
// create the rotations
for (int j = 0; j < (COLS - 1); j++) {
double a = matrix[j][j];
double b = matrix[j + 1][j];
cblas_drotg(&a, &b, &c[j], &s[j]);
}
// apply the rotations along the diagonal
char side = 'L';
char pivot = 'V';
char direct = 'F';
lapack_int m = ROWS;
lapack_int n = COLS;
lapack_int lda = ROWS;
dlasr(&side, &pivot, &direct, &m, &n, c, s, (double*) matrix, &lda);
print_matrix();
return 0;
}