1 /*******************************************************************************
2 * PRIMME PReconditioned Iterative MultiMethod Eigensolver
3 * Copyright (C) 2018 College of William & Mary,
4 * James R. McCombs, Eloy Romero Alcalde, Andreas Stathopoulos, Lingfei Wu
5 *
6 * This file is part of PRIMME.
7 *
8 * PRIMME is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 *
13 * PRIMME is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with this library; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
21 *
22 *******************************************************************************
23 * File: rsbw.c
24 *
25 * Purpose - librsb wrapper.
26 *
27 ******************************************************************************/
28
29 #include <rsb.h> /* for rsb_lib_init */
30 #include <blas_sparse.h>
31 #include <stdio.h>
32 #include <assert.h>
33 #include "rsbw.h"
34 #include "num.h"
35 #include <math.h>
36
readMatrixRSB(const char * matrixFileName,blas_sparse_matrix * matrix,double * fnorm)37 int readMatrixRSB(const char* matrixFileName, blas_sparse_matrix *matrix, double *fnorm) {
38 #if !defined(RSB_NUMERICAL_TYPE_DOUBLE) || !defined(RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX)
39 fprintf(stderr, "Needed librsb with support for 'double' and 'double complex'.\n");
40 return -1;
41 #else
42 # ifdef USE_DOUBLECOMPLEX
43 rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX;
44 # else
45 rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE;
46 # endif
47 if (rsb_lib_init(RSB_NULL_INIT_OPTIONS) != RSB_ERR_NO_ERROR) assert(0);
48 *matrix = blas_invalid_handle;
49 if ((rsb_perror(NULL, rsb_lib_init(RSB_NULL_INIT_OPTIONS)))!=RSB_ERR_NO_ERROR) {
50 fprintf(stderr, "Error while initializing librsb.\n");
51 return -1;
52 }
53
54 *matrix = rsb_load_spblas_matrix_file_as_matrix_market(matrixFileName, typecode);
55 if ( *matrix == blas_invalid_handle) {
56 fprintf(stderr, "ERROR: Could not read matrix file\n");
57 return -1;
58 }
59
60 if (BLAS_ussp(*matrix, blas_rsb_autotune_next_operation) != 0) assert(0);
61 if (BLAS_dusget_infinity_norm(*matrix, fnorm, blas_no_trans) != 0) assert(0);
62
63 return 0;
64 #endif
65 }
66
RSBMatvec(void * x,void * y,int * blockSize,primme_params * primme)67 void RSBMatvec(void *x, void *y, int *blockSize, primme_params *primme) {
68 int i;
69 SCALAR *xvec, *yvec;
70 blas_sparse_matrix *matrix;
71 #ifdef USE_DOUBLECOMPLEX
72 const SCALAR one=(SCALAR)1;
73 #endif
74
75 if (*blockSize <= 0) return;
76
77 matrix = (blas_sparse_matrix *)primme->matrix;
78 xvec = (SCALAR *)x;
79 yvec = (SCALAR *)y;
80
81 for (i=0; i<*blockSize*primme->nLocal; i++)
82 yvec[i] = 0;
83 #ifndef USE_DOUBLECOMPLEX
84 if(BLAS_dusmm(blas_colmajor, blas_no_trans, *blockSize, 1.0, *matrix, xvec, primme->nLocal, yvec, primme->nLocal) != 0) assert(0);
85 #else
86 if(BLAS_zusmm(blas_colmajor, blas_no_trans, *blockSize, &one, *matrix, xvec, primme->nLocal, yvec, primme->nLocal) != 0) assert(0);
87 #endif
88 }
89
RSBMatvecSVD(void * x,int * ldx,void * y,int * ldy,int * blockSize,int * trans,primme_svds_params * primme_svds)90 void RSBMatvecSVD(void *x, int *ldx, void *y, int *ldy, int *blockSize,
91 int *trans, primme_svds_params *primme_svds) {
92
93 int i, j;
94 SCALAR *xvec, *yvec;
95 blas_sparse_matrix *matrix;
96 #ifdef USE_DOUBLECOMPLEX
97 const SCALAR one=(SCALAR)1;
98 #endif
99
100 matrix = (blas_sparse_matrix *)primme_svds->matrix;
101 xvec = (SCALAR *)x;
102 yvec = (SCALAR *)y;
103
104 for (i=0; i<(*blockSize); i++) {
105 if (*trans == 0){
106 for (j=0; j<primme_svds->mLocal; j++) yvec[i*(*ldy)+j] = 0;
107 } else {
108 for (j=0; j<primme_svds->nLocal; j++) yvec[i*(*ldy)+j] = 0;
109 }
110 }
111
112 #ifndef USE_DOUBLECOMPLEX
113 if (BLAS_dusmm(blas_colmajor, *trans == 0 ? blas_no_trans : blas_trans, *blockSize, 1.0, *matrix, xvec, *ldx, yvec, *ldy) != 0) assert(0);
114 #else
115 if (BLAS_zusmm(blas_colmajor, *trans == 0 ? blas_no_trans : blas_trans, *blockSize, &one, *matrix, xvec, *ldx, yvec, *ldy) != 0) assert(0);
116 #endif
117 }
118
119
120 /******************************************************************************
121 * Generates the diagonal of A.
122 *
123 * P = Diag(A)
124 *
125 * This will be used with solver provided shifts as (P-shift_i)^(-1)
126 ******************************************************************************/
getDiagonal(blas_sparse_matrix matrix,double * diag)127 static void getDiagonal(blas_sparse_matrix matrix, double *diag) {
128 #ifndef USE_DOUBLECOMPLEX
129 if (BLAS_dusget_diag(matrix, diag) != 0) assert(0);
130 #else
131 int n = BLAS_usgp(matrix, blas_num_rows), i;
132 SCALAR *d = (SCALAR *)primme_calloc(n, sizeof(PRIMM_NUM), "aux");
133 if (BLAS_zusget_diag(matrix, d) != 0) assert(0);
134 for (i=0; i<n; i++) diag[i] = REAL_PART(d[i]);
135 #endif
136 }
137
createInvDiagPrecRSB(blas_sparse_matrix matrix,double shift,double ** prec)138 int createInvDiagPrecRSB(blas_sparse_matrix matrix, double shift, double **prec) {
139 int i, n = BLAS_usgp(matrix, blas_num_rows);
140 double *diag;
141
142 diag = (double*)primme_calloc(n, sizeof(double), "diag");
143 getDiagonal(matrix, diag);
144 for (i=0; i<n; i++)
145 diag[i] -= shift;
146 *prec = diag;
147 return 1;
148 }
149
150
151 /******************************************************************************
152 * Generates sum of square values per rows and then per columns
153 *
154 ******************************************************************************/
getSumSquares(blas_sparse_matrix matrix,double * diag)155 static void getSumSquares(blas_sparse_matrix matrix, double *diag) {
156 int nnz = BLAS_usgp(matrix, blas_num_nonzeros);
157 int m = BLAS_usgp(matrix, blas_num_rows);
158 int n = BLAS_usgp(matrix, blas_num_cols);
159 int i, *AI, *AJ;
160 double *sumr = diag, *sumc = &diag[m], v;
161 SCALAR *A;
162
163 for (i=0; i < m + n; i++) {
164 diag[i] = 0.0;
165 }
166
167 A = (SCALAR *)primme_calloc(nnz, sizeof(SCALAR), "A");
168 AI = (int *)primme_calloc(nnz*2, sizeof(int), "AI AJ");
169 AJ = AI + nnz;
170 for (i=0; i<nnz; i++) {
171 v = A[i]*CONJ(A[i]);
172 sumr[AI[i]] += v;
173 sumc[AJ[i]] += v;
174 }
175 }
176
createInvNormalPrecRSB(blas_sparse_matrix matrix,double shift,double ** prec)177 int createInvNormalPrecRSB(blas_sparse_matrix matrix, double shift, double **prec) {
178 int i;
179 double *diag, minDenominator=1e-14;
180 int n = BLAS_usgp(matrix, blas_num_rows)+BLAS_usgp(matrix, blas_num_cols);
181
182 diag = (double*)primme_calloc(n, sizeof(double), "diag");
183 getSumSquares(matrix, diag);
184 for (i=0; i<n; i++) {
185 diag[i] -= shift*shift;
186 if (fabs(diag[i]) < minDenominator)
187 diag[i] = copysign(minDenominator, diag[i]);
188 }
189 *prec = diag;
190 return 1;
191 }
192