power.c File Reference
Example programs and code

A toy <rsb.h>-based C program implementing the power method for computing matrix eigenvalues. Uses rsb_spmv(). More...

## Functions

int | main (const int argc, char *const argv[])

## Detailed Description

A toy <rsb.h>-based C program implementing the power method for computing matrix eigenvalues. Uses rsb_spmv().

Author
  Michele Martone

/*
Copyright (C) 2008-2021 Michele Martone
This file is part of librsb.
librsb is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
librsb 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 Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with librsb; see the file COPYING.
If not, see <https://www.gnu.org/licenses/>.
*/
/*!
@file
@author Michele Martone
@brief A toy <rsb.h>-based C program implementing the power
method for computing matrix eigenvalues. Uses #rsb_spmv().
\ingroup rsb_doc_examples
\include power.c
*/
#include <stdio.h> // printf
#include <math.h> // sqrt
#include <stdlib.h> // calloc
#include <rsb.h>
int main(const int argc, char * const argv[])
{
const int WANT_VERBOSE = 0;
struct rsb_mtx_t *mtxAp = NULL;
const int bs = RSB_DEFAULT_BLOCKING;
const int br = bs, bc = bs; /* bs x bs blocked */
const rsb_nnz_idx_t nnzA = 4;
const rsb_coo_idx_t nrA = 3;
const rsb_coo_idx_t ncA = 3;
rsb_int_t it = 0;
const rsb_int_t maxit = 100;
const rsb_coo_idx_t IA[] = { 0, 1, 2, 0 };
const rsb_coo_idx_t JA[] = { 0, 1, 2, 2 };
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE VA[] = { 11, 22, 33, 13 };
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE ZERO = 0;
int i;
rsb_err_t errval = 0;
RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE norm = 0.0, /* nu */
oldnorm = 1.0, /* oldnorm */
*b1 = NULL, *b2 = NULL,
*bnow = NULL, *bnext = NULL;/* b1 and b2 aliases */
rsb_type_t typecode = RSB_NUMERICAL_TYPE_FIRST_BLAS;
size_t ds = 0;
/* tolerance */
const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE tol = 1e-14;
/* library initialization */
if(rsb_lib_init(RSB_NULL_INIT_OPTIONS)!=RSB_ERR_NO_ERROR)
return EXIT_FAILURE;
/* allocation */
mtxAp = rsb_mtx_alloc_from_coo_const(VA,IA,JA,nnzA,
typecode,nrA,ncA,br,bc,RSB_FLAG_NOFLAGS,NULL);
if(!mtxAp)
return EXIT_FAILURE;
ds = (nrA)*sizeof(RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE);
b1 = calloc(1,ds);
b2 = calloc(1,ds);
if(! (b1 && b2))
{
errval = RSB_ERR_ENOMEM;
goto err;
}
for( i = 0; i < nrA; ++i )
b1[i] = 1;
bnow = b1, bnext = b2;/* b,b' */
while( fabs(norm-oldnorm) > tol && it<maxit )
{
++ it;
oldnorm = norm;
/* b'<-Ab */
if(( rsb_spmv(RSB_TRANSPOSITION_N,NULL,mtxAp,bnow,
1,&ZERO,bnext,1)) != RSB_ERR_NO_ERROR )
goto err;
/* nu<-||Ab||^2 */
norm = 0;
for(i=0;i<nrA;++i)
norm += bnext[i]*bnext[i];
/* nu<-||Ab|| */
norm = sqrt(norm);
norm = 1.0/norm;
/* b'<- Ab / ||Ab|| */
for(i=0;i<nrA;++i)
bnext[i] *= norm;
norm = 1.0/norm;
printf("it:%d norm:%lg norm diff:%lg\n",it,norm,norm-oldnorm);
{void *tmp=bnow;bnow=bnext;bnext=tmp;/* pointers swap */}
if(WANT_VERBOSE)
{
printf("norm:%lg\n",norm);
if(isinf(norm))
/* isinf is a C99 feature (need correct
* compilation flags) */
goto err;
for(i=0;i<2;++i)
printf("x[%d]=%lg\n",i,((double*)bnext)[i]);
}
}
/* the biggest eigenvalue should be in bnow */
rsb_mtx_free(mtxAp);
free(b1);
free(b2);
if(rsb_lib_exit(RSB_NULL_EXIT_OPTIONS)!=RSB_ERR_NO_ERROR)
goto err;
if( it == maxit )
{
printf("ERROR: hit iterations limit without convergence!");
errval=RSB_ERR_GENERIC_ERROR;
}
return EXIT_SUCCESS;
err:
rsb_perror(NULL,errval);
return EXIT_FAILURE;
}
main
auto main() -> int
Definition assemble.cpp:38
rsb_spmv
rsb_err_t rsb_spmv(rsb_trans_t transA, const void *alphap, const struct rsb_mtx_t *mtxAp, const void *Xp, rsb_coo_idx_t incX, const void *betap, void *Yp, rsb_coo_idx_t incY)
Definition rsb_rsb.c:529
rsb_lib_exit
rsb_err_t rsb_lib_exit(struct rsb_initopts *iop)
Definition rsb_rsb.c:202
rsb_perror
rsb_err_t rsb_perror(void *stream, rsb_err_t errval)
Definition rsb_rsb.c:814
rsb_mtx_free
struct rsb_mtx_t * rsb_mtx_free(struct rsb_mtx_t *mtxAp)
Definition rsb_rsb.c:301
rsb_lib_init
rsb_err_t rsb_lib_init(struct rsb_initopts *iop)
Definition rsb_rsb.c:56
rsb_mtx_alloc_from_coo_const
struct rsb_mtx_t * rsb_mtx_alloc_from_coo_const(const void *VA, const rsb_coo_idx_t *IA, const rsb_coo_idx_t *JA, rsb_nnz_idx_t nnzA, rsb_type_t typecode, rsb_coo_idx_t nrA, rsb_coo_idx_t ncA, rsb_blk_idx_t brA, rsb_blk_idx_t bcA, rsb_flags_t flagsA, rsb_err_t *errvalp)
Definition rsb_rsb.c:240
rsb.h
This file declares the user interface functions and data structures for the librsb library (see The l...
RSB_FLAG_NOFLAGS
#define RSB_FLAG_NOFLAGS
Definition rsb.h:478
RSB_ERR_NO_ERROR
#define RSB_ERR_NO_ERROR
Definition rsb.h:616
RSB_NULL_EXIT_OPTIONS
#define RSB_NULL_EXIT_OPTIONS
Definition rsb.h:848
RSB_DEFAULT_BLOCKING
#define RSB_DEFAULT_BLOCKING
Definition rsb.h:429
rsb_nnz_idx_t
signed int rsb_nnz_idx_t
Definition rsb.h:362
rsb_coo_idx_t
signed int rsb_coo_idx_t
Definition rsb.h:349
RSB_ERR_ENOMEM
#define RSB_ERR_ENOMEM
Definition rsb.h:637
rsb_err_t
signed int rsb_err_t
Definition rsb.h:386
rsb_type_t
char rsb_type_t
Definition rsb.h:378
RSB_ERR_GENERIC_ERROR
#define RSB_ERR_GENERIC_ERROR
Definition rsb.h:619
RSB_NULL_INIT_OPTIONS
#define RSB_NULL_INIT_OPTIONS
Definition rsb.h:847
rsb_int_t
signed int rsb_int_t
Definition rsb.h:392
RSB_TRANSPOSITION_N
#define RSB_TRANSPOSITION_N
N: Non transposed flag, valid for rsb_trans_t typed variables.
Definition rsb_types.h:138
RSB_NUMERICAL_TYPE_FIRST_BLAS
#define RSB_NUMERICAL_TYPE_FIRST_BLAS
A default numerical matrix type; if possible, not integer one. If no such type is configured in,...
Definition rsb_types.h:203
RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE
#define RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE
The default, blas if possible , numerical type (can be used for declarations).
Definition rsb_types.h:75

## Function Documentation

## main()

int main | ( | const int | argc,
|  | char *const | argv[] )