diff --git a/include/spm.h b/include/spm.h
index fbe39e71e925c31540c3a598feaf2393334e6393..58a9689dd1a10cc6289702d0319a7a4ef20e7145 100644
--- a/include/spm.h
+++ b/include/spm.h
@@ -48,9 +48,9 @@
*/
typedef struct spmatrix_s {
spm_mtxtype_t mtxtype; /**< Matrix structure: SpmGeneral, SpmSymmetric
- or SpmHermitian. */
+ or SpmHermitian. */
spm_coeftype_t flttype; /**< values datatype: SpmPattern, SpmFloat, SpmDouble,
- SpmComplex32 or SpmComplex64 */
+ SpmComplex32 or SpmComplex64 */
spm_fmttype_t fmttype; /**< Matrix storage format: SpmCSC, SpmCSR, SpmIJV */
spm_int_t gN; /**< Global number of vertices in the compressed graph (Computed) */
@@ -64,8 +64,8 @@ typedef struct spmatrix_s {
spm_int_t nnzexp; /**< Local number of non zeroes in the compressed graph (Computed) */
spm_int_t dof; /**< Number of degrees of freedom per unknown,
- if > 0, constant degree of freedom
- otherwise, irregular degree of freedom (refer to dofs) */
+ if > 0, constant degree of freedom
+ otherwise, irregular degree of freedom (refer to dofs) */
spm_int_t *dofs; /**< Array of the first column of each element in the
expanded matrix [+baseval] */
spm_layout_t layout; /**< SpmColMajor, or SpmRowMajor */
@@ -98,9 +98,9 @@ void spmGenFakeValues( spmatrix_t *spm );
double spmNorm( spm_normtype_t ntype, const spmatrix_t *spm );
int spmMatVec( spm_trans_t trans, const void *alpha, const spmatrix_t *spm, const void *x, const void *beta, void *y );
int spmMatMat( spm_trans_t trans, spm_int_t n,
- const void *alpha, const spmatrix_t *A,
- const void *B, spm_int_t ldb,
- const void *beta, void *C, spm_int_t ldc );
+ const void *alpha, const spmatrix_t *A,
+ const void *B, spm_int_t ldb,
+ const void *beta, void *C, spm_int_t ldc );
void spmScalMatrix( double alpha, spmatrix_t *spm );
void spmScalVector( spm_coeftype_t flt, double alpha, spm_int_t n, void *x, spm_int_t incx );
@@ -146,9 +146,9 @@ int spmSave( const spmatrix_t *spm, FILE *outfile );
* @{
*/
int spmReadDriver( spm_driver_t driver,
- const char *filename,
- spmatrix_t *spm,
- MPI_Comm spm_comm );
+ const char *filename,
+ spmatrix_t *spm,
+ MPI_Comm spm_comm );
/**
* @}
* @name SPM debug subroutines
diff --git a/src/integer_sort_mtypes.c b/src/integer_sort_mtypes.c
index b1662a00881c1b5aefc1d0c674f1d15d1736ffd4..882974fa05cb4e3990c763df762be2fe00a173e0 100644
--- a/src/integer_sort_mtypes.c
+++ b/src/integer_sort_mtypes.c
@@ -78,10 +78,10 @@ typedef struct
log(MAX_THRESH_2)). Since total_elements has type size_t, we get as
upper bound for log (total_elements):
bits per unsigned char (CHAR_BIT) * sizeof(size_t). */
-#define STACK_SIZE_2 (CHAR_BIT * sizeof (spm_int_t))
-#define PUSH_2(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
-#define POP_2(low, high) ((void) (--top, (low = top->lo), (high = top->hi)))
-#define STACK_NOT_EMPTY_2 (stack < top)
+#define STACK_SIZE_2 (CHAR_BIT * sizeof (spm_int_t))
+#define PUSH_2(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
+#define POP_2(low, high) ((void) (--top, (low = top->lo), (high = top->hi)))
+#define STACK_NOT_EMPTY_2 (stack < top)
#endif /* MAX_THRESH_2 */
diff --git a/src/z_spm.c b/src/z_spm.c
index 2dc33b8ab7d6ad9565da62ff4c9f86863a7cb7b3..37f534d7a45c9329161efba68ca2d5523a49d1e8 100644
--- a/src/z_spm.c
+++ b/src/z_spm.c
@@ -50,42 +50,42 @@ z_spmSort( spmatrix_t *spm )
(void)sortptr;
if (spm->dof > 1){
- fprintf(stderr, "z_spmSort: Number of dof (%d) greater than one not supported\n", (int)spm->dof);
- exit(1);
+ fprintf(stderr, "z_spmSort: Number of dof (%d) greater than one not supported\n", (int)spm->dof);
+ exit(1);
}
/* Sort in place each subset */
if ( spm->fmttype == SpmCSC ) {
- for (i=0; i<n; i++, colptr++)
- {
- size = colptr[1] - colptr[0];
+ for (i=0; i<n; i++, colptr++)
+ {
+ size = colptr[1] - colptr[0];
#if defined(PRECISION_p)
- spmIntSort1Asc1( rowptr, size );
+ spmIntSort1Asc1( rowptr, size );
#else
- sortptr[0] = rowptr;
- sortptr[1] = values;
- z_spmIntFltSortAsc( sortptr, size );
+ sortptr[0] = rowptr;
+ sortptr[1] = values;
+ z_spmIntFltSortAsc( sortptr, size );
#endif
- rowptr += size;
- values += size;
- }
+ rowptr += size;
+ values += size;
+ }
}
else if ( spm->fmttype == SpmCSR ) {
- for (i=0; i<n; i++, rowptr++)
- {
- size = rowptr[1] - rowptr[0];
+ for (i=0; i<n; i++, rowptr++)
+ {
+ size = rowptr[1] - rowptr[0];
#if defined(PRECISION_p)
- spmIntSort1Asc1( rowptr, size );
+ spmIntSort1Asc1( rowptr, size );
#else
- sortptr[0] = colptr;
- sortptr[1] = values;
- z_spmIntFltSortAsc( sortptr, size );
+ sortptr[0] = colptr;
+ sortptr[1] = values;
+ z_spmIntFltSortAsc( sortptr, size );
#endif
- colptr += size;
- values += size;
- }
+ colptr += size;
+ values += size;
+ }
}
}
@@ -128,55 +128,55 @@ z_spmMergeDuplicate( spmatrix_t *spm )
spm_int_t merge = 0;
if ( spm->fmttype == SpmCSC ) {
- idx = 0;
- for (i=0; i<n; i++, colptr++)
- {
- size = colptr[1] - colptr[0];
- for (j = 0; j < size;
- j++, oldrow++, oldval+=dof2, newrow++, newval+=dof2, idx++)
- {
- if ( newrow != oldrow ) {
- newrow[0] = oldrow[0];
+ idx = 0;
+ for (i=0; i<n; i++, colptr++)
+ {
+ size = colptr[1] - colptr[0];
+ for (j = 0; j < size;
+ j++, oldrow++, oldval+=dof2, newrow++, newval+=dof2, idx++)
+ {
+ if ( newrow != oldrow ) {
+ newrow[0] = oldrow[0];
#if !defined(PRECISION_p)
- memcpy( newval, oldval, dof2 * sizeof(spm_complex64_t) );
+ memcpy( newval, oldval, dof2 * sizeof(spm_complex64_t) );
#endif
- }
+ }
- while( ((j+1) < size) && (newrow[0] == oldrow[1]) ) {
- j++; oldrow++; oldval+=dof2;
+ while( ((j+1) < size) && (newrow[0] == oldrow[1]) ) {
+ j++; oldrow++; oldval+=dof2;
#if !defined(PRECISION_p)
- /* Merge the two sets of values */
- for (d = 0; d < dof2; d++) {
- newval[d] += oldval[d];
- }
+ /* Merge the two sets of values */
+ for (d = 0; d < dof2; d++) {
+ newval[d] += oldval[d];
+ }
#endif
- merge++;
- }
- }
- assert( ((merge == 0) && ( colptr[1] == idx+baseval)) ||
- ((merge != 0) && ( colptr[1] > idx+baseval)) );
-
- colptr[1] = idx + baseval;
- }
- assert( ((merge == 0) && (spm->nnz == idx)) ||
- ((merge != 0) && (spm->nnz - merge == idx)) );
-
- if (merge > 0) {
- spm->nnz = spm->nnz - merge;
- spm->gnnz = spm->nnz;
-
- newrow = malloc( spm->nnz * sizeof(spm_int_t) );
- memcpy( newrow, spm->rowptr, spm->nnz * sizeof(spm_int_t) );
- free(spm->rowptr);
- spm->rowptr = newrow;
+ merge++;
+ }
+ }
+ assert( ((merge == 0) && ( colptr[1] == idx+baseval)) ||
+ ((merge != 0) && ( colptr[1] > idx+baseval)) );
+
+ colptr[1] = idx + baseval;
+ }
+ assert( ((merge == 0) && (spm->nnz == idx)) ||
+ ((merge != 0) && (spm->nnz - merge == idx)) );
+
+ if (merge > 0) {
+ spm->nnz = spm->nnz - merge;
+ spm->gnnz = spm->nnz;
+
+ newrow = malloc( spm->nnz * sizeof(spm_int_t) );
+ memcpy( newrow, spm->rowptr, spm->nnz * sizeof(spm_int_t) );
+ free(spm->rowptr);
+ spm->rowptr = newrow;
#if !defined(PRECISION_p)
- newval = malloc( spm->nnz * dof2 * sizeof(spm_int_t) );
- memcpy( newval, spm->values, spm->nnz * dof2 * sizeof(spm_complex64_t) );
- free(spm->values);
- spm->values = newval;
+ newval = malloc( spm->nnz * dof2 * sizeof(spm_int_t) );
+ memcpy( newval, spm->values, spm->nnz * dof2 * sizeof(spm_complex64_t) );
+ free(spm->values);
+ spm->values = newval;
#endif
- }
+ }
}
(void)d;
@@ -223,153 +223,153 @@ z_spmSymmetrize( spmatrix_t *spm )
toaddsze = 0;
if ( (spm->fmttype == SpmCSC) || (spm->fmttype == SpmCSR) ) {
- if (spm->fmttype == SpmCSC) {
- oldcol = spm->colptr;
- coltmp = spm->colptr;
- oldrow = spm->rowptr;
- rowtmp = spm->rowptr;
- }
- else {
- oldcol = spm->rowptr;
- coltmp = spm->rowptr;
- oldrow = spm->colptr;
- rowtmp = spm->colptr;
- }
-
- baseval = oldcol[0];
- for (i=0; i<n; i++, coltmp++)
- {
- size = coltmp[1] - coltmp[0];
- for (r=0; r<size; r++, rowtmp++ )
- {
- j = rowtmp[0]-baseval;
- if ( i != j ) {
- /* Look for the element (j, i) */
- spm_int_t frow = oldcol[ j ] - baseval;
- spm_int_t lrow = oldcol[ j+1 ] - baseval;
- int found = 0;
-
- for (k = frow; (k < lrow); k++)
- {
- if (i == (oldrow[k]-baseval))
- {
- found = 1;
- break;
- }
- else if ( i < (oldrow[k]-baseval))
- {
- /* The spm is sorted so we will not find it later */
- break;
- }
- }
-
- if ( !found ) {
- if ( newcol == NULL ) {
- newcol = malloc( (spm->n+1) * sizeof(spm_int_t) );
- for (k=0; k<spm->n; k++) {
- newcol[k] = oldcol[k+1] - oldcol[k];
- }
-
- /* Let's start with a buffer that will contain 5% of extra elements */
- toaddsze = spm_imax( 0.05 * (double)spm->nnz, 1 );
- toaddtab = malloc( 2*toaddsze * sizeof(spm_int_t) );
- }
-
- if (toaddcnt >= toaddsze) {
- toaddsze *= 2;
- toaddtab = (spm_int_t*)realloc( toaddtab, 2*toaddsze*sizeof(spm_int_t) );
- }
-
- /* Newcol stores the number of elements per column */
- newcol[ j ]++;
- /* toaddtab stores the couples (j, i) to be added in the final sparse matrix */
- toaddtab[ toaddcnt * 2 ] = j;
- toaddtab[ toaddcnt * 2 + 1 ] = i;
- toaddcnt++;
- }
- }
- }
- }
-
- if (toaddcnt > 0) {
- spm_int_t newsize, oldsize;
-
- /* Sort the array per column */
- spmIntSort2Asc1(toaddtab, toaddcnt);
-
- spm->nnz = spm->nnz + toaddcnt;
- spm->gnnz = spm->nnz;
-
- newrow = malloc( spm->nnz * sizeof(spm_int_t) );
+ if (spm->fmttype == SpmCSC) {
+ oldcol = spm->colptr;
+ coltmp = spm->colptr;
+ oldrow = spm->rowptr;
+ rowtmp = spm->rowptr;
+ }
+ else {
+ oldcol = spm->rowptr;
+ coltmp = spm->rowptr;
+ oldrow = spm->colptr;
+ rowtmp = spm->colptr;
+ }
+
+ baseval = oldcol[0];
+ for (i=0; i<n; i++, coltmp++)
+ {
+ size = coltmp[1] - coltmp[0];
+ for (r=0; r<size; r++, rowtmp++ )
+ {
+ j = rowtmp[0]-baseval;
+ if ( i != j ) {
+ /* Look for the element (j, i) */
+ spm_int_t frow = oldcol[ j ] - baseval;
+ spm_int_t lrow = oldcol[ j+1 ] - baseval;
+ int found = 0;
+
+ for (k = frow; (k < lrow); k++)
+ {
+ if (i == (oldrow[k]-baseval))
+ {
+ found = 1;
+ break;
+ }
+ else if ( i < (oldrow[k]-baseval))
+ {
+ /* The spm is sorted so we will not find it later */
+ break;
+ }
+ }
+
+ if ( !found ) {
+ if ( newcol == NULL ) {
+ newcol = malloc( (spm->n+1) * sizeof(spm_int_t) );
+ for (k=0; k<spm->n; k++) {
+ newcol[k] = oldcol[k+1] - oldcol[k];
+ }
+
+ /* Let's start with a buffer that will contain 5% of extra elements */
+ toaddsze = spm_imax( 0.05 * (double)spm->nnz, 1 );
+ toaddtab = malloc( 2*toaddsze * sizeof(spm_int_t) );
+ }
+
+ if (toaddcnt >= toaddsze) {
+ toaddsze *= 2;
+ toaddtab = (spm_int_t*)realloc( toaddtab, 2*toaddsze*sizeof(spm_int_t) );
+ }
+
+ /* Newcol stores the number of elements per column */
+ newcol[ j ]++;
+ /* toaddtab stores the couples (j, i) to be added in the final sparse matrix */
+ toaddtab[ toaddcnt * 2 ] = j;
+ toaddtab[ toaddcnt * 2 + 1 ] = i;
+ toaddcnt++;
+ }
+ }
+ }
+ }
+
+ if (toaddcnt > 0) {
+ spm_int_t newsize, oldsize;
+
+ /* Sort the array per column */
+ spmIntSort2Asc1(toaddtab, toaddcnt);
+
+ spm->nnz = spm->nnz + toaddcnt;
+ spm->gnnz = spm->nnz;
+
+ newrow = malloc( spm->nnz * sizeof(spm_int_t) );
#if !defined(PRECISION_p)
- newval = malloc( spm->nnz * dof2 * sizeof(spm_complex64_t) );
+ newval = malloc( spm->nnz * dof2 * sizeof(spm_complex64_t) );
#endif
- coltmp = newcol;
- rowtmp = newrow;
- valtmp = newval;
- oldval = spm->values;
- toaddtmp = toaddtab;
-
- newsize = coltmp[0];
- coltmp[0] = baseval;
- for (i=0; i<n; i++, coltmp++, oldcol++)
- {
- /* Copy the existing elements */
- oldsize = oldcol[1] - oldcol[0];
- memcpy( rowtmp, oldrow, oldsize * sizeof(spm_int_t) );
+ coltmp = newcol;
+ rowtmp = newrow;
+ valtmp = newval;
+ oldval = spm->values;
+ toaddtmp = toaddtab;
+
+ newsize = coltmp[0];
+ coltmp[0] = baseval;
+ for (i=0; i<n; i++, coltmp++, oldcol++)
+ {
+ /* Copy the existing elements */
+ oldsize = oldcol[1] - oldcol[0];
+ memcpy( rowtmp, oldrow, oldsize * sizeof(spm_int_t) );
#if !defined(PRECISION_p)
- memcpy( valtmp, oldval, oldsize * dof2 * sizeof(spm_complex64_t) );
+ memcpy( valtmp, oldval, oldsize * dof2 * sizeof(spm_complex64_t) );
#endif
- oldrow += oldsize;
- rowtmp += oldsize;
- oldval += oldsize * dof2;
- valtmp += oldsize * dof2;
+ oldrow += oldsize;
+ rowtmp += oldsize;
+ oldval += oldsize * dof2;
+ valtmp += oldsize * dof2;
- /* Some elements have been added */
- assert( newsize >= oldsize );
- if ( newsize > oldsize ) {
- int added = 0;
- int tobeadded = newsize - oldsize;
+ /* Some elements have been added */
+ assert( newsize >= oldsize );
+ if ( newsize > oldsize ) {
+ int added = 0;
+ int tobeadded = newsize - oldsize;
- /* At least one element is equal to i */
- assert( toaddtmp[0] == i );
+ /* At least one element is equal to i */
+ assert( toaddtmp[0] == i );
- /* Let's add the new vertices */
- while( (added < tobeadded) && (toaddtmp[0] == i) ) {
- rowtmp[0] = toaddtmp[1] + baseval;
- rowtmp++; toaddtmp+=2; added++;
- }
- assert( added == tobeadded );
+ /* Let's add the new vertices */
+ while( (added < tobeadded) && (toaddtmp[0] == i) ) {
+ rowtmp[0] = toaddtmp[1] + baseval;
+ rowtmp++; toaddtmp+=2; added++;
+ }
+ assert( added == tobeadded );
#if !defined(PRECISION_p)
- /* Add 0.s for the associated values */
- memset( valtmp, 0, added * dof2 * sizeof(spm_complex64_t) );
- valtmp += added * dof2;
+ /* Add 0.s for the associated values */
+ memset( valtmp, 0, added * dof2 * sizeof(spm_complex64_t) );
+ valtmp += added * dof2;
#endif
- }
-
- /* Use oldsize as temporary variable to update the new colptr */
- oldsize = newsize;
- newsize = coltmp[1];
- coltmp[1] = coltmp[0] + oldsize;
- }
-
- assert( coltmp[0]-baseval == spm->nnz );
- free( toaddtab );
- free( spm->colptr );
- free( spm->rowptr );
- free( spm->values );
- if (spm->fmttype == SpmCSC) {
- spm->colptr = newcol;
- spm->rowptr = newrow;
- }
- else {
- spm->colptr = newrow;
- spm->rowptr = newcol;
- }
- spm->values = newval;
- }
+ }
+
+ /* Use oldsize as temporary variable to update the new colptr */
+ oldsize = newsize;
+ newsize = coltmp[1];
+ coltmp[1] = coltmp[0] + oldsize;
+ }
+
+ assert( coltmp[0]-baseval == spm->nnz );
+ free( toaddtab );
+ free( spm->colptr );
+ free( spm->rowptr );
+ free( spm->values );
+ if (spm->fmttype == SpmCSC) {
+ spm->colptr = newcol;
+ spm->rowptr = newrow;
+ }
+ else {
+ spm->colptr = newrow;
+ spm->rowptr = newcol;
+ }
+ spm->values = newval;
+ }
}
return toaddcnt;
diff --git a/src/z_spm_integer.c b/src/z_spm_integer.c
index f0ef07bddbbbe3517418c761538a43a9b81799bb..760a999a30ff63f4c41aad4b143d3f1e6584498a 100644
--- a/src/z_spm_integer.c
+++ b/src/z_spm_integer.c
@@ -94,24 +94,24 @@ static size_t intsortsize_iif[3] = { sizeof(spm_int_t), sizeof(spm_int_t), sizeo
#define INTSORTSIZE(x) (intsortsize_iif[x])
#define INTSORTNTAB 3
#define INTSORTSWAP(p,q) do { \
- spm_int_t t; \
- long disp_p = (((spm_int_t*)p)-((spm_int_t*)base_ptr)); \
- long disp_q = (((spm_int_t*)q)-((spm_int_t*)base_ptr)); \
- spm_int_t * intptr = *(pbase+1); \
- spm_complex64_t * fltptr = *(pbase+2); \
- spm_complex64_t f; \
- /* swap integers */ \
- t = *((spm_int_t *) (p)); \
- *((spm_int_t *) (p)) = *((spm_int_t *) (q)); \
- *((spm_int_t *) (q)) = t; \
- /* swap on second integer array */ \
- t = intptr[disp_p]; \
- intptr[disp_p] = intptr[disp_q]; \
- intptr[disp_q] = t; \
- /* swap corresponding values */ \
- f = fltptr[disp_p]; \
- fltptr[disp_p] = fltptr[disp_q]; \
- fltptr[disp_q] = f; \
+ spm_int_t t; \
+ long disp_p = (((spm_int_t*)p)-((spm_int_t*)base_ptr)); \
+ long disp_q = (((spm_int_t*)q)-((spm_int_t*)base_ptr)); \
+ spm_int_t * intptr = *(pbase+1); \
+ spm_complex64_t * fltptr = *(pbase+2); \
+ spm_complex64_t f; \
+ /* swap integers */ \
+ t = *((spm_int_t *) (p)); \
+ *((spm_int_t *) (p)) = *((spm_int_t *) (q)); \
+ *((spm_int_t *) (q)) = t; \
+ /* swap on second integer array */ \
+ t = intptr[disp_p]; \
+ intptr[disp_p] = intptr[disp_q]; \
+ intptr[disp_q] = t; \
+ /* swap corresponding values */ \
+ f = fltptr[disp_p]; \
+ fltptr[disp_p] = fltptr[disp_q]; \
+ fltptr[disp_q] = f; \
} while (0)
static inline int
diff --git a/src/z_spm_matrixvector.c b/src/z_spm_matrixvector.c
index 9396cd45efe55520b50469cfa88f72216c9f5498..f0fc6c63dd8076a456ad0afab425820952a887aa 100644
--- a/src/z_spm_matrixvector.c
+++ b/src/z_spm_matrixvector.c
@@ -63,11 +63,11 @@
*******************************************************************************/
int
z_spmGeCSCv(const spm_trans_t trans,
- spm_complex64_t alpha,
- const spmatrix_t *spm,
- const spm_complex64_t *x,
- spm_complex64_t beta,
- spm_complex64_t *y )
+ spm_complex64_t alpha,
+ const spmatrix_t *spm,
+ const spm_complex64_t *x,
+ spm_complex64_t beta,
+ spm_complex64_t *y )
{
const spm_complex64_t *valptr = (spm_complex64_t*)(spm->values);
const spm_complex64_t *xptr = (const spm_complex64_t*)x;
@@ -76,73 +76,73 @@ z_spmGeCSCv(const spm_trans_t trans,
if ( (spm == NULL) || (x == NULL) || (y == NULL ) )
{
- return SPM_ERR_BADPARAMETER;
+ return SPM_ERR_BADPARAMETER;
}
if( spm->mtxtype != SpmGeneral )
{
- return SPM_ERR_BADPARAMETER;
+ return SPM_ERR_BADPARAMETER;
}
baseval = spmFindBase( spm );
/* first, y = beta*y */
if( beta == 0. ) {
- memset( yptr, 0, spm->gN * sizeof(spm_complex64_t) );
+ memset( yptr, 0, spm->gN * sizeof(spm_complex64_t) );
}
else {
- for( i=0; i<spm->gN; i++, yptr++ ) {
- (*yptr) *= beta;
- }
- yptr = y;
+ for( i=0; i<spm->gN; i++, yptr++ ) {
+ (*yptr) *= beta;
+ }
+ yptr = y;
}
if( alpha != 0. ) {
- /**
- * SpmNoTrans
- */
- if( trans == SpmNoTrans )
- {
- for( col=0; col < spm->gN; col++ )
- {
- for( i=spm->colptr[col]; i<spm->colptr[col+1]; i++ )
- {
- row = spm->rowptr[i-baseval]-baseval;
- yptr[row] += alpha * valptr[i-baseval] * xptr[col];
- }
- }
- }
- /**
- * SpmTrans
- */
- else if( trans == SpmTrans )
- {
- for( col=0; col < spm->gN; col++ )
- {
- for( i=spm->colptr[col]; i<spm->colptr[col+1]; i++ )
- {
- row = spm->rowptr[i-baseval]-baseval;
- yptr[col] += alpha * valptr[i-baseval] * xptr[row];
- }
- }
- }
+ /**
+ * SpmNoTrans
+ */
+ if( trans == SpmNoTrans )
+ {
+ for( col=0; col < spm->gN; col++ )
+ {
+ for( i=spm->colptr[col]; i<spm->colptr[col+1]; i++ )
+ {
+ row = spm->rowptr[i-baseval]-baseval;
+ yptr[row] += alpha * valptr[i-baseval] * xptr[col];
+ }
+ }
+ }
+ /**
+ * SpmTrans
+ */
+ else if( trans == SpmTrans )
+ {
+ for( col=0; col < spm->gN; col++ )
+ {
+ for( i=spm->colptr[col]; i<spm->colptr[col+1]; i++ )
+ {
+ row = spm->rowptr[i-baseval]-baseval;
+ yptr[col] += alpha * valptr[i-baseval] * xptr[row];
+ }
+ }
+ }
#if defined(PRECISION_c) || defined(PRECISION_z)
- else if( trans == SpmConjTrans )
- {
- for( col=0; col < spm->gN; col++ )
- {
- for( i=spm->colptr[col]; i<spm->colptr[col+1]; i++ )
- {
- row = spm->rowptr[i-baseval]-baseval;
- yptr[col] += alpha * conj( valptr[i-baseval] ) * xptr[row];
- }
- }
- }
+ else if( trans == SpmConjTrans )
+ {
+ for( col=0; col < spm->gN; col++ )
+ {
+ for( i=spm->colptr[col]; i<spm->colptr[col+1]; i++ )
+ {
+ row = spm->rowptr[i-baseval]-baseval;
+ yptr[col] += alpha * conj( valptr[i-baseval] ) * xptr[row];
+ }
+ }
+ }
#endif
- else
- {
- return SPM_ERR_BADPARAMETER;
- }
+ else
+ {
+ return SPM_ERR_BADPARAMETER;
+ }
}
return SPM_SUCCESS;
@@ -184,10 +184,10 @@ z_spmGeCSCv(const spm_trans_t trans,
*******************************************************************************/
int
z_spmSyCSCv( spm_complex64_t alpha,
- const spmatrix_t *spm,
- const spm_complex64_t *x,
- spm_complex64_t beta,
- spm_complex64_t *y )
+ const spmatrix_t *spm,
+ const spm_complex64_t *x,
+ spm_complex64_t beta,
+ spm_complex64_t *y )
{
const spm_complex64_t *valptr = (spm_complex64_t*)spm->values;
const spm_complex64_t *xptr = x;
@@ -196,40 +196,40 @@ z_spmSyCSCv( spm_complex64_t alpha,
if ( (spm == NULL) || (x == NULL) || (y == NULL ) )
{
- return SPM_ERR_BADPARAMETER;
+ return SPM_ERR_BADPARAMETER;
}
if( spm->mtxtype != SpmSymmetric )
{
- return SPM_ERR_BADPARAMETER;
+ return SPM_ERR_BADPARAMETER;
}
baseval = spmFindBase( spm );
/* First, y = beta*y */
if( beta == 0. ) {
- memset( yptr, 0, spm->gN * sizeof(spm_complex64_t) );
+ memset( yptr, 0, spm->gN * sizeof(spm_complex64_t) );
}
else {
- for( i=0; i<spm->gN; i++, yptr++ ) {
- (*yptr) *= beta;
- }
- yptr = y;
+ for( i=0; i<spm->gN; i++, yptr++ ) {
+ (*yptr) *= beta;
+ }
+ yptr = y;
}
if( alpha != 0. ) {
- for( col=0; col < spm->gN; col++ )
- {
- for( i=spm->colptr[col]; i < spm->colptr[col+1]; i++ )
- {
- row = spm->rowptr[i-baseval]-baseval;
- yptr[row] += alpha * valptr[i-baseval] * xptr[col];
- if( col != row )
- {
- yptr[col] += alpha * valptr[i-baseval] * xptr[row];
- }
- }
- }
+ for( col=0; col < spm->gN; col++ )
+ {
+ for( i=spm->colptr[col]; i < spm->colptr[col+1]; i++ )
+ {
+ row = spm->rowptr[i-baseval]-baseval;
+ yptr[row] += alpha * valptr[i-baseval] * xptr[col];
+ if( col != row )
+ {
+ yptr[col] += alpha * valptr[i-baseval] * xptr[row];
+ }
+ }
+ }
}
return SPM_SUCCESS;
@@ -272,10 +272,10 @@ z_spmSyCSCv( spm_complex64_t alpha,
*******************************************************************************/
int
z_spmHeCSCv( spm_complex64_t alpha,
- const spmatrix_t *spm,
- const spm_complex64_t *x,
- spm_complex64_t beta,
- spm_complex64_t *y )
+ const spmatrix_t *spm,
+ const spm_complex64_t *x,
+ spm_complex64_t beta,
+ spm_complex64_t *y )
{
const spm_complex64_t *valptr = (spm_complex64_t*)spm->values;
const spm_complex64_t *xptr = x;
@@ -284,42 +284,42 @@ z_spmHeCSCv( spm_complex64_t alpha,
if ( (spm == NULL) || (x == NULL) || (y == NULL ) )
{
- return SPM_ERR_BADPARAMETER;
+ return SPM_ERR_BADPARAMETER;
}
if( spm->mtxtype != SpmHermitian )
{
- return SPM_ERR_BADPARAMETER;
+ return SPM_ERR_BADPARAMETER;
}
/* First, y = beta*y */
if( beta == 0. ) {
- memset( yptr, 0, spm->gN * sizeof(spm_complex64_t) );
+ memset( yptr, 0, spm->gN * sizeof(spm_complex64_t) );
}
else {
- for( i=0; i<spm->gN; i++, yptr++ ) {
- (*yptr) *= beta;
- }
- yptr = y;
+ for( i=0; i<spm->gN; i++, yptr++ ) {
+ (*yptr) *= beta;
+ }
+ yptr = y;
}
baseval = spmFindBase( spm );
if( alpha != 0. ) {
- for( col=0; col < spm->gN; col++ )
- {
- for( i=spm->colptr[col]; i < spm->colptr[col+1]; i++ )
- {
- row=spm->rowptr[i-baseval]-baseval;
- if( col != row ) {
- yptr[row] += alpha * valptr[i-baseval] * xptr[col];
- yptr[col] += alpha * conj( valptr[i-baseval] ) * xptr[row];
- }
- else {
- yptr[row] += alpha * creal(valptr[i-baseval]) * xptr[col];
- }
- }
- }
+ for( col=0; col < spm->gN; col++ )
+ {
+ for( i=spm->colptr[col]; i < spm->colptr[col+1]; i++ )
+ {
+ row=spm->rowptr[i-baseval]-baseval;
+ if( col != row ) {
+ yptr[row] += alpha * valptr[i-baseval] * xptr[col];
+ yptr[col] += alpha * conj( valptr[i-baseval] ) * xptr[row];
+ }
+ else {
+ yptr[row] += alpha * creal(valptr[i-baseval]) * xptr[col];
+ }
+ }
+ }
}
return SPM_SUCCESS;
@@ -365,11 +365,11 @@ z_spmHeCSCv( spm_complex64_t alpha,
*******************************************************************************/
int
z_spmCSCMatVec(const spm_trans_t trans,
- const void *alphaptr,
- const spmatrix_t *spm,
- const void *xptr,
- const void *betaptr,
- void *yptr )
+ const void *alphaptr,
+ const spmatrix_t *spm,
+ const void *xptr,
+ const void *betaptr,
+ void *yptr )
{
const spm_complex64_t *x = (const spm_complex64_t*)xptr;
spm_complex64_t *y = (spm_complex64_t*)yptr;
@@ -381,13 +381,13 @@ z_spmCSCMatVec(const spm_trans_t trans,
switch (spm->mtxtype) {
#if defined(PRECISION_z) || defined(PRECISION_c)
case SpmHermitian:
- return z_spmHeCSCv( alpha, spm, x, beta, y );
+ return z_spmHeCSCv( alpha, spm, x, beta, y );
#endif
case SpmSymmetric:
- return z_spmSyCSCv( alpha, spm, x, beta, y );
+ return z_spmSyCSCv( alpha, spm, x, beta, y );
case SpmGeneral:
default:
- return z_spmGeCSCv( trans, alpha, spm, x, beta, y );
+ return z_spmGeCSCv( trans, alpha, spm, x, beta, y );
}
}
@@ -446,14 +446,14 @@ z_spmCSCMatVec(const spm_trans_t trans,
*******************************************************************************/
int
z_spmCSCMatMat(const spm_trans_t trans,
- spm_int_t n,
- const void *alphaptr,
- const spmatrix_t *A,
- const void *Bptr,
- spm_int_t ldb,
- const void *betaptr,
- void *Cptr,
- spm_int_t ldc )
+ spm_int_t n,
+ const void *alphaptr,
+ const spmatrix_t *A,
+ const void *Bptr,
+ spm_int_t ldb,
+ const void *betaptr,
+ void *Cptr,
+ spm_int_t ldc )
{
const spm_complex64_t *B = (const spm_complex64_t*)Bptr;
spm_complex64_t *C = (spm_complex64_t*)Cptr;
@@ -466,21 +466,21 @@ z_spmCSCMatMat(const spm_trans_t trans,
switch (A->mtxtype) {
#if defined(PRECISION_z) || defined(PRECISION_c)
case SpmHermitian:
- for( i=0; i<n; i++ ){
- rc = z_spmHeCSCv( alpha, A, B + i * ldb, beta, C + i *ldc );
- }
- break;
+ for( i=0; i<n; i++ ){
+ rc = z_spmHeCSCv( alpha, A, B + i * ldb, beta, C + i *ldc );
+ }
+ break;
#endif
case SpmSymmetric:
- for( i=0; i<n; i++ ){
- rc = z_spmSyCSCv( alpha, A, B + i * ldb, beta, C + i *ldc );
- }
- break;
+ for( i=0; i<n; i++ ){
+ rc = z_spmSyCSCv( alpha, A, B + i * ldb, beta, C + i *ldc );
+ }
+ break;
case SpmGeneral:
default:
- for( i=0; i<n; i++ ){
- rc = z_spmGeCSCv( trans, alpha, A, B + i * ldb, beta, C + i *ldc );
- }
+ for( i=0; i<n; i++ ){
+ rc = z_spmGeCSCv( trans, alpha, A, B + i * ldb, beta, C + i *ldc );
+ }
}
return rc;
}