1 /* pdnmesh - a 2D finite element solver
2 Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>
3 This program is free software; you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation; either version 2 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 $Id: eig_arpack.c,v 1.7 2005/04/21 23:16:54 sarod Exp $
17 */
18 /* ARPACK Solver */
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22
23 #include <math.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26 #ifdef HAVE_PTHREAD
27 #include <pthread.h>
28 #endif
29 #include "types.h"
30
31 int arpack_max_iterations=300;
32 int arpack_subspace_dimension=20;
33
34 /* prototype for the fortran function */
35 extern void ardrive_(int *N, int *Nz, double *AA, int *JA, int *IA,
36 int *nev, int *ncv, int *lworkl, double *resid,
37 double *v, double *workd, double *workl, int *info, int *ido,
38 double *select_, double *d, int *ierr,
39 int *NL, double *LL, int *JL, int *IL,
40 int *NU, double *UU, int *JU, int *IU, int *maxit);
41
42 #ifdef HAVE_PTHREAD
43 typedef struct t_barrier_ {
44 int tcount; /* current no. of threads inside barrier */
45 int nthreads; /* the no. of threads the barrier works
46 with. This is a constant */
47 pthread_mutex_t enter_mutex;
48 pthread_mutex_t exit_mutex;
49 pthread_cond_t lastthread_cond;
50 pthread_cond_t exit_cond;
51 } th_barrier;
52
53 /* initialize barrier */
54 /* N - no. of accomodated threads */
55 static void
init_th_barrier(th_barrier * barrier,int N)56 init_th_barrier(th_barrier *barrier, int N)
57 {
58 barrier->tcount=0; /* initially empty */
59 barrier->nthreads=N;
60 pthread_mutex_init(&barrier->enter_mutex,NULL);
61 pthread_mutex_init(&barrier->exit_mutex,NULL);
62 pthread_cond_init(&barrier->lastthread_cond,NULL);
63 pthread_cond_init(&barrier->exit_cond,NULL);
64 }
65 /* destroy barrier */
66 static void
destroy_th_barrier(th_barrier * barrier)67 destroy_th_barrier(th_barrier *barrier)
68 {
69 pthread_mutex_destroy(&barrier->enter_mutex);
70 pthread_mutex_destroy(&barrier->exit_mutex);
71 pthread_cond_destroy(&barrier->lastthread_cond);
72 pthread_cond_destroy(&barrier->exit_cond);
73 barrier->tcount=barrier->nthreads=0;
74 }
75 /* the main operation of the barrier */
76 static void
sync_barrier(th_barrier * barrier)77 sync_barrier(th_barrier *barrier)
78 {
79 /* trivial case */
80 if(barrier->nthreads <2) return;
81 /* else */
82 /* new threads enters the barrier. Now close the entry door
83 so that other threads cannot enter the barrier until we are done */
84 pthread_mutex_lock(&barrier->enter_mutex);
85 /* next lock the exit mutex - no threads can leave either */
86 pthread_mutex_lock(&barrier->exit_mutex);
87 /* now check to see if this is the last expected thread */
88 if( ++(barrier->tcount) < barrier->nthreads) {
89 /* no. this is not the last thread. so open the entry door */
90 pthread_mutex_unlock(&barrier->enter_mutex);
91 /* go to sleep */
92 pthread_cond_wait(&barrier->exit_cond,&barrier->exit_mutex);
93 } else {
94 /* this is the last thread */
95 /* wakeup sleeping threads */
96 pthread_cond_broadcast(&barrier->exit_cond);
97 /* go to sleep until all threads are woken up
98 and leave the barrier */
99 pthread_cond_wait(&barrier->lastthread_cond,&barrier->exit_mutex);
100 /* now all threads have left the barrier. so open the entry door again */
101 pthread_mutex_unlock(&barrier->enter_mutex);
102 }
103 /* next to the last thread leaving the barrier */
104 if(--(barrier->tcount)==1) {
105 /* wakeup the last sleeping thread */
106 pthread_cond_broadcast(&barrier->lastthread_cond);
107 }
108 pthread_mutex_unlock(&barrier->exit_mutex);
109 }
110
111 /* structure to pass to the LU threads */
112 typedef struct tdata_ {
113 th_barrier *barrier;
114 int myid;
115 int *N;
116 SPMat *A;
117 SPMat *L;
118 int *T; /* no. of threads */
119 pthread_mutex_t *lock_mutex; /* lock for updating L*/
120 } tdata;
121
122 static void *
work_threadfn(void * data)123 work_threadfn(void *data)
124 {
125 tdata *t=(tdata*)data;
126 int c; /* current column */
127 int i,k;
128 double temp,ll,aa,bb;
129 int N=SPM_dim(t->A);
130
131 for (c=0; c<N; c++) {
132 #ifdef DEBUG
133 printf("col %d, enter id=%d\n",c,t->myid);
134 #endif
135 /* the 0-th thread will calculate the diagonal entry */
136 if(t->myid==0) {
137 temp=0;
138 /* diagonal element of row j */
139 for (i=0;i<c; i++) {
140 aa=SPM_e(t->L,c,i);
141 if(aa)
142 temp+=aa*aa;
143 }
144 temp=SPM_e(t->A,c,c)-temp;
145
146 if(temp<=0) {
147 fprintf(stderr,"%s: %d: matrix not positive definite\n",__FILE__,__LINE__);
148 temp=TOL;
149 } else {
150 temp=sqrt(temp);
151 }
152 pthread_mutex_lock(t->lock_mutex);
153 SPM_eq(t->L,c,c,temp);
154 pthread_mutex_unlock(t->lock_mutex);
155 #ifdef DEBUG
156 printf("thread %d, L[%d][%d]=%lf\n",t->myid,c,c,temp);
157 #endif
158 }
159
160 sync_barrier(t->barrier);
161 /* other threads will calculate the remaining entries of column c */
162 if(t->myid>0) {
163 i=c+t->myid;
164 ll=SPM_e(t->L,c,c);
165 while (i<N) {
166 temp=0;
167 for(k=0; k<c; k++) {
168 aa=SPM_e(t->L,i,k);
169 bb=SPM_e(t->L,c,k);
170 if((aa)&&(bb))
171 temp+=aa*bb;
172 #ifdef DEBUG
173 printf("[%d][%d] add [%d][%d]*[%d][%d]\n",i,c,i,k,c,k);
174 #endif
175 }
176 aa=(SPM_e(t->A,i,c)-temp)/ll;
177 pthread_mutex_lock(t->lock_mutex);
178 SPM_eq(t->L,i,c,aa);
179 pthread_mutex_unlock(t->lock_mutex);
180 #ifdef DEBUG
181 printf("thread %d, L[%d][%d]=%lf\n",t->myid,i,c,aa);
182 #endif
183 i+=(*(t->T)-1);
184 }
185 }
186 sync_barrier(t->barrier);
187 #ifdef DEBUG
188 printf("col %d, leave id=%d\n",c,t->myid);
189 #endif
190 }
191 return NULL;
192 }
193
194 /* Cholevsky decpmposition with thread support */
195 static int
cholevsky_decompose(SPMat * A,SPMat * L)196 cholevsky_decompose(SPMat *A, SPMat *L) {
197 int M=3;
198 int i;
199 tdata *thread_data;
200 th_barrier b;
201 pthread_t *th_array;
202 pthread_attr_t attr;
203 pthread_mutex_t lock_mutex;
204
205 /* threads =3, one master, two slave */
206 /* thread related things */
207 init_th_barrier(&b,M);
208 pthread_attr_init(&attr);
209 pthread_attr_setdetachstate(&attr,PTHREAD_CREATE_JOINABLE);
210 pthread_mutex_init(&lock_mutex,NULL);
211
212
213 if ((th_array=(pthread_t*)malloc((size_t)M*sizeof(pthread_t)))==0) {
214 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
215 exit(1);
216 }
217 if ((thread_data=(tdata*)malloc((size_t)M*sizeof(tdata)))==0) {
218 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
219 exit(1);
220 }
221
222 for (i=0; i<M; i++) {
223 thread_data[i].myid=i;
224 thread_data[i].barrier=&b;
225 thread_data[i].T=&M;
226 thread_data[i].A=A;
227 thread_data[i].L=L;
228 thread_data[i].lock_mutex=&lock_mutex;
229 pthread_create(&th_array[i],&attr,work_threadfn,(void*)(&thread_data[i]));
230 }
231
232 for(i=0; i<M; i++) {
233 pthread_join(th_array[i],NULL);
234 }
235
236 destroy_th_barrier(&b);
237 pthread_attr_destroy(&attr);
238 pthread_mutex_destroy(&lock_mutex);
239 free(th_array);
240 free(thread_data);
241
242 return 0;
243 }
244 #endif /* ! HAVE_PTHREAD */
245
246 #ifndef HAVE_PTHREAD
247 /* create cholevsky decomposition of sparse matrix A,
248 output L, the sparse matrix.
249 L should be initialized first
250 */
251 static int
cholevsky_decompose(SPMat * A,SPMat * L)252 cholevsky_decompose(SPMat *A, SPMat *L) {
253 int i,j,k,N;
254 double aa,bb,ll,temp;
255 N=SPM_dim(A);
256 if(SPM_dim(L)!=N) return -1;
257
258 for(j=0; j<N; j++) {
259 temp=0;
260 /* diagonal element of row j */
261 for (k=0;k<j; k++) {
262 aa=SPM_e(L,j,k);
263 if(aa)
264 temp+=aa*aa;
265 }
266 temp=SPM_e(A,j,j)-temp;
267
268 if(temp<=0) {
269 fprintf(stderr,"%s: %d: matrix not positive definite\n",__FILE__,__LINE__);
270 temp=TOL;
271 } else {
272 temp=sqrt(temp);
273 }
274 SPM_eq(L,j,j,temp);
275 ll=temp; /* remember this value */
276 /* now find the j th column below the diagonal element */
277 for (i=j+1; i<N; i++) {
278 temp=0;
279 for(k=0; k<j; k++) {
280 aa=SPM_e(L,i,k);
281 bb=SPM_e(L,j,k);
282 if((aa)&&(bb))
283 temp+=aa*bb;
284 }
285 aa=(SPM_e(A,i,j)-temp)/ll;
286 SPM_eq(L,i,j,aa);
287 }
288 }
289
290 return 0;
291 }
292 #endif /* HAVE_PTHREAD */
293
294 /* AA - CSR format lower triangular matrix */
295 /* solve Ax = y using forward elimination */
296 static int
forward_eliminate(double * AA,int * JA,int * IA,int N,int Nz,double * x,double * y)297 forward_eliminate(double *AA, int *JA, int *IA, int N, int Nz,
298 double *x, double *y) {
299 double temp,ll;
300 int i,k,k1,k2;
301 for(i=0; i<N; i++) {
302 k1=IA[i];
303 if (i<N-1) {
304 k2=IA[i+1]-1;
305 } else {
306 k2=Nz-1;
307 }
308 /* range of row i is [k1,k2] */
309 /* find the product AA(k1:k2) (except diagonal) with x(0:i-1) */
310 /* note we have no ordering of row elements,
311 so we do not know where the diagonal element is */
312 ll=1;
313 temp=0;
314 for(k=k1; k<=k2; k++) {
315 if(JA[k]!=i) { /* avoid diagonal element */
316 temp+=AA[k]*x[JA[k]];
317 } else {
318 ll=AA[k];
319 }
320 }
321 x[i]=(y[i]-temp)/ll;
322 }
323
324 return 0;
325 }
326 /* solve Ax = y by back substitution, A is an upper triangular
327 matrix, given in CSR format */
328 static int
backward_substitute(double * AA,int * JA,int * IA,int N,int Nz,double * x,double * y)329 backward_substitute(double *AA, int *JA, int *IA, int N, int Nz,
330 double *x, double *y) {
331 double temp,ll;
332 int i,k,k1,k2;
333
334 for (i=N-1; i>=0; i--) {
335 k1=IA[i];
336 if (i<N-1) {
337 k2=IA[i+1]-1;
338 } else {
339 k2=Nz-1;
340 }
341 /* range for row i is [k1,k2]
342 find the product of row i, AA(k1:k2) except diagonal
343 with x[i+1:N-1] */
344 temp=0;
345 ll=1;
346 for(k=k1;k<=k2; k++) {
347 if(JA[k]!=i) {
348 temp+=AA[k]*x[JA[k]];
349 } else {
350 ll=AA[k];
351 }
352 }
353 x[i]=(y[i]-temp)/ll;
354 }
355 return 0;
356 }
357
358 /* find generalized eigenvalues of the pencil (A,B) using ARPACK
359 store eigenvectors in E - size N x nev
360 and eigenvalues in r - size nev x 1
361 nev - the no. of eigenvalues to find
362 lower_tiangular - if this is 1,
363 NOTE: A, B are lower triangular, so need to fill up the upper triangle
364 of the matrices. Otherwise, full matrices are assumed.
365 */
366 int
arpack_find_eigs(SPMat * A,SPMat * B,MY_DOUBLE ** E,MY_DOUBLE * r,MY_INT nev,int lower_triangular)367 arpack_find_eigs(SPMat *A, SPMat *B, MY_DOUBLE **E, MY_DOUBLE *r, MY_INT nev, int lower_triangular) {
368 SPMat L,U;
369 int N,i,j;
370
371 double *AA,*AL, *AU;
372 int *IA, *JA, *IL, *JL, *IU, *JU;
373 unsigned int M,LL,UU;
374 int ncv, lworkl;
375 double *workd, *workl, *v, *resid;
376 double *selection, *d;
377 int info, ido, ierr;
378
379 double *x, *y, temp;
380 int maxiter;
381
382 N=SPM_dim(A);
383
384 if( lower_triangular ) {
385 /* fill up the upper triangle of both matrices, first */
386 for (i=0; i<N; i++) {
387 for (j=0; j<i; j++) {
388 temp=SPM_e(A,i,j);
389 if(temp)
390 SPM_eq(A,j,i,temp);
391 temp=SPM_e(B,i,j);
392 if(temp)
393 SPM_eq(B,j,i,temp);
394 }
395 }
396 }
397
398 /* find cholevsky decomposition of B*/
399 SPM_init(&L,N);
400 cholevsky_decompose(B, &L);
401 #ifdef DEBUG
402 printf("A=[\n");
403 for (i=0; i<N; i++) {
404 for(j=0; j<N; j++) {
405 printf("%lf ",SPM_e(A,i,j));
406 }
407 printf("\n");
408 }
409 printf("];\n");
410 #endif
411 #ifdef DEBUG
412 printf("B=[\n");
413 for (i=0; i<N; i++) {
414 for(j=0; j<N; j++) {
415 printf("%lf ",SPM_e(B,i,j));
416 }
417 printf("\n");
418 }
419 printf("];\n");
420 #endif
421
422
423 #ifdef DEBUG
424 printf("L=[\n");
425 for (i=0; i<N; i++) {
426 for(j=0; j<N; j++) {
427 printf("%lf ",SPM_e(&L,i,j));
428 }
429 printf("\n");
430 }
431 printf("];\n");
432 #endif
433 printf("A: sparsity=%lf, B: sparsity=%lf, L: sparsity=%lf, size=%d\n",
434 (double)SPM_size(A)/(double)(N*N),(double)SPM_size(B)/(double)(N*N),
435 (double)SPM_size(&L)/(double)(N*N),N);
436 M=SPM_size(A);
437 /* CSR arrays of A */
438 if ((AA=(double *)calloc((size_t)M,sizeof(double)))==0) {
439 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
440 exit(1);
441 }
442 if ((JA=(int *)calloc((size_t)M,sizeof(int)))==0) {
443 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
444 exit(1);
445 }
446 if ((IA=(int *)calloc((size_t)N,sizeof(int)))==0) {
447 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
448 exit(1);
449 }
450
451 SPM_build_lists(A);
452 info=SPM_build_CSR(A,AA,JA,IA);
453 #ifdef DEBUG
454 printf("A (%d) size=%d, CSR=%d\n",N,M,info);
455 #endif
456 /* convert JA and IA to fortran format */
457 for(i=0; i<N; i++) {
458 JA[i]+=1;
459 IA[i]+=1;
460 }
461 for(i=N; i<M; i++)
462 JA[i]+=1;
463 #ifdef DEBUG1
464 for(i=0; i<M; i++) {
465 if(i<N) {
466 printf("%lf : %d : %d\n",AA[i],JA[i],IA[i]);
467 } else {
468 printf("%lf : %d\n",AA[i],JA[i]);
469 }
470 }
471 #endif
472
473 LL=SPM_size(&L);
474 /* CSR arrays */
475 if ((AL=(double *)calloc((size_t)LL,sizeof(double)))==0) {
476 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
477 exit(1);
478 }
479 if ((JL=(int *)calloc((size_t)LL,sizeof(int)))==0) {
480 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
481 exit(1);
482 }
483 if ((IL=(int *)calloc((size_t)N,sizeof(int)))==0) {
484 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
485 exit(1);
486 }
487
488 SPM_build_lists(&L);
489 info=SPM_build_CSR(&L,AL,JL,IL);
490 #ifdef DEBUG1
491 printf("L size=%d, CSR=%d\n",LL,info);
492 for(i=0; i<LL; i++) {
493 if(i<N) {
494 printf("%lf : %d : %d\n",AL[i],JL[i],IL[i]);
495 } else {
496 printf("%lf : %d\n",AL[i],JL[i]);
497 }
498 }
499 #endif
500
501 /* create L^T as U */
502 SPM_init(&U,N);
503 SPM_transpose(&L,&U);
504 SPM_destroy(&L);
505
506 SPM_build_lists(&U);
507 UU=SPM_size(&U);
508 /* CSR arrays */
509 if ((AU=(double *)calloc((size_t)UU,sizeof(double)))==0) {
510 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
511 exit(1);
512 }
513 if ((JU=(int *)calloc((size_t)UU,sizeof(int)))==0) {
514 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
515 exit(1);
516 }
517 if ((IU=(int *)calloc((size_t)N,sizeof(int)))==0) {
518 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
519 exit(1);
520 }
521
522 info=SPM_build_CSR(&U,AU,JU,IU);
523
524 #ifdef DEBUG
525 printf("L^T size=%d, CSR=%d\n",UU,info);
526 #endif
527 SPM_destroy(&U);
528 #ifdef DEBUG1
529 for(i=0; i<UU; i++) {
530 if(i<N) {
531 printf("%lf : %d : %d\n",AU[i],JU[i],IU[i]);
532 } else {
533 printf("%lf : %d\n",AU[i],JU[i]);
534 }
535 }
536 #endif
537
538 /* convert arrays to fortran format */
539 for(i=0; i<N; i++) {
540 JU[i]+=1;
541 JL[i]+=1;
542 IU[i]+=1;
543 IL[i]+=1;
544 }
545 for(i=N; i<UU; i++) {
546 JU[i]+=1;
547 JL[i]+=1;
548 }
549
550 maxiter=arpack_max_iterations;
551 if(maxiter<0) maxiter=300;
552
553 /* typically we set the subspace dimension to 4*nev,
554 However, if we want a higher value, we use the global variable */
555 ncv=4*nev;
556 if(ncv <arpack_subspace_dimension) {
557 ncv=arpack_subspace_dimension;
558 } else {
559 arpack_subspace_dimension=ncv;
560 }
561 if(ncv<4) ncv=4;
562 if(ncv>N) ncv=N-1;
563
564 lworkl=ncv*(ncv+8);
565 /* ARPACK storage allocation */
566 if ((workd=(double *)calloc((size_t)3*N,sizeof(double)))==0) {
567 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
568 exit(1);
569 }
570 if ((workl=(double *)calloc((size_t)lworkl,sizeof(double)))==0) {
571 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
572 exit(1);
573 }
574 if ((v=(double *)calloc((size_t)N*ncv,sizeof(double)))==0) {
575 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
576 exit(1);
577 }
578 if ((resid=(double *)calloc((size_t)N,sizeof(double)))==0) {
579 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
580 exit(1);
581 }
582 if ((selection=(double *)calloc((size_t)ncv,sizeof(double)))==0) {
583 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
584 exit(1);
585 }
586 if ((d=(double *)calloc((size_t)ncv*2,sizeof(double)))==0) {
587 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
588 exit(1);
589 }
590
591 ardrive_(&N, &M, AA, JA, IA,
592 &nev, &ncv, &lworkl, resid,
593 v, workd, workl, &info, &ido,
594 selection, d, &ierr,
595 &LL, AL, JL, IL, &UU, AU, JU, IU, &maxiter);
596
597 printf("info=%d, ido=%d, ierr=%d\n",info,ido,ierr);
598
599 free(AA);
600 free(JA);
601 free(IA);
602 free(AL);
603 free(JL);
604 free(IL);
605 free(workd);
606 free(workl);
607 free(resid);
608 free(selection);
609
610 /* convert arrays back to C format */
611 for(i=0; i<N; i++) {
612 JU[i]-=1;
613 IU[i]-=1;
614 }
615 for(i=N; i<UU; i++) {
616 JU[i]-=1;
617 }
618 if ((x=(double *)calloc((size_t)N,sizeof(double)))==0) {
619 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
620 exit(1);
621 }
622 if ((y=(double *)calloc((size_t)N,sizeof(double)))==0) {
623 fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
624 exit(1);
625 }
626 for(i=0; i<nev; i++){
627 r[i]=d[i];
628 for(j=0; j<N; j++)
629 y[j]=v[i*N+j];
630 backward_substitute(AU, JU, IU, N, UU, x, y);
631 for(j=0; j<N; j++)
632 E[j][i]=x[j];
633 }
634
635 free(d);
636 free(v);
637 free(AU);
638 free(IU);
639 free(JU);
640
641
642 free(x);
643 free(y);
644 return 0;
645 }
646
647 /* dummy drivers when we do not use ARPACK */
648 #ifndef USE_ARPACK
dsaupd_(int * ido,char * bmat,int * n,char * which,int * nev,double * tol,double * resid,int * ncv,double * v,int * ldv,int * param,int * ipntr,double * workd,double * workl,int * lworkl,int * info,unsigned long bmat_len,unsigned long which_len)649 void dsaupd_(int *ido, char *bmat, int *n, char *which,
650 int *nev, double *tol, double *resid, int *ncv,
651 double *v, int *ldv, int *param, int *ipntr,
652 double *workd, double *workl, int *lworkl,
653 int *info, unsigned long bmat_len, unsigned long which_len) {
654 *ido=0;
655 *info=0;
656 }
657
dseupd_(int * rvec,char * all,int * select_,double * d,double * v1,int * ldv1,double * sigma,char * bmat,int * n,char * which,int * nev,double * tol,double * resid,int * ncv,double * v,int * ldv,int * iparam,int * ipntr,double * workd,double * workl,int * lworkl,int * ierr)658 void dseupd_(int *rvec, char *all, int *select_, double *d,
659 double *v1, int *ldv1, double *sigma,
660 char *bmat, int *n, char *which, int *nev,
661 double *tol, double *resid, int *ncv, double *v,
662 int *ldv, int *iparam, int *ipntr, double *workd,
663 double *workl, int *lworkl, int *ierr) {
664 *ierr=0;
665 }
666 #endif /* USE_ARPACK */
667