Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
PreAMGInterpEM.c
Go to the documentation of this file.
1
20#include <math.h>
21#include <time.h>
22
23#ifdef _OPENMP
24#include <omp.h>
25#endif
26
27#include "fasp.h"
28#include "fasp_functs.h"
29
30/*---------------------------------*/
31/*-- Declare Private Functions --*/
32/*---------------------------------*/
33
34static SHORT getiteval(dCSRmat *, dCSRmat *);
35static SHORT invden(INT, REAL *, REAL *);
36static SHORT get_block(dCSRmat *, INT, INT, INT *, INT *, REAL *, INT *);
37static SHORT gentisquare_nomass(dCSRmat *, INT, INT *, REAL *, INT *);
38static SHORT getinonefull(INT **, REAL **, INT *, INT, INT *, REAL *);
39static SHORT orderone(INT **, REAL **, INT *);
40static SHORT genintval(dCSRmat *, INT **, REAL **, INT, INT *, INT, INT, INT);
41
42/*---------------------------------*/
43/*-- Public Functions --*/
44/*---------------------------------*/
45
63void fasp_amg_interp_em (dCSRmat *A,
64 ivector *vertices,
65 dCSRmat *P,
66 AMG_param *param)
67{
68 INT *vec = vertices->val;
69 INT *CoarseIndex = (INT *)fasp_mem_calloc(vertices->row, sizeof(INT));
70 INT i, j, index;
71
72 // generate indices for C-points
73 for ( index = i = 0; i < vertices->row; ++i ) {
74 if ( vec[i] == 1 ) {
75 CoarseIndex[i] = index;
76 index++;
77 }
78 }
79
80#ifdef _OPENMP
81#pragma omp parallel for private(i,j) if(P->nnz>OPENMP_HOLDS)
82#endif
83 for ( i = 0; i < P->nnz; ++i ) {
84 j = P->JA[i];
85 P->JA[i] = CoarseIndex[j];
86 }
87
88 // clean up memory
89 fasp_mem_free(CoarseIndex); CoarseIndex = NULL;
90
91 // main part
92 getiteval(A, P);
93}
94
95/*---------------------------------*/
96/*-- Private Functions --*/
97/*---------------------------------*/
98
115static SHORT invden (INT nn,
116 REAL *mat,
117 REAL *invmat)
118{
119 INT i,j;
120 SHORT status = FASP_SUCCESS;
121
122#ifdef _OPENMP
123 // variables for OpenMP
124 INT myid, mybegin, myend;
125 INT nthreads = fasp_get_num_threads();
126#endif
127
128 INT *pivot=(INT *)fasp_mem_calloc(nn,sizeof(INT));
129 REAL *rhs=(REAL *)fasp_mem_calloc(nn,sizeof(REAL));
130 REAL *sol=(REAL *)fasp_mem_calloc(nn,sizeof(REAL));
131
132 fasp_smat_lu_decomp(mat,pivot,nn);
133
134#ifdef _OPENMP
135#pragma omp parallel for private(myid,mybegin,myend,i,j) if(nn>OPENMP_HOLDS)
136 for (myid=0; myid<nthreads; ++myid) {
137 fasp_get_start_end(myid, nthreads, nn, &mybegin, &myend);
138 for (i=mybegin; i<myend; ++i) {
139#else
140 for (i=0;i<nn;++i) {
141#endif
142 for (j=0;j<nn;++j) rhs[j]=0.;
143 rhs[i]=1.;
144 fasp_smat_lu_solve(mat,rhs,pivot,sol,nn);
145 for (j=0;j<nn;++j) invmat[i*nn+j]=sol[j];
146#ifdef _OPENMP
147 }
148 }
149#else
150 }
151#endif
152
153 fasp_mem_free(pivot); pivot = NULL;
154 fasp_mem_free(rhs); rhs = NULL;
155 fasp_mem_free(sol); sol = NULL;
156
157 return status;
158}
159
181static SHORT get_block (dCSRmat *A,
182 INT m,
183 INT n,
184 INT *rows,
185 INT *cols,
186 REAL *Aloc,
187 INT *mask)
188{
189 INT i, j, k, iloc;
190
191#ifdef _OPENMP
192 // variables for OpenMP
193 INT myid, mybegin, myend;
194 INT nthreads = fasp_get_num_threads();
195#endif
196
197 memset(Aloc, 0x0, sizeof(REAL)*m*n);
198
199#ifdef _OPENMP
200#pragma omp parallel for if(n>OPENMP_HOLDS) private(j)
201#endif
202 for ( j=0; j<n; ++j ) {
203 mask[cols[j]] = j; // initialize mask, mask stores C indices 0,1,...
204 }
205
206#ifdef _OPENMP
207#pragma omp parallel for private(myid,mybegin,myend,i,j,k,iloc) if(m>OPENMP_HOLDS)
208 for ( myid=0; myid<nthreads; ++myid ) {
209 fasp_get_start_end(myid, nthreads, m, &mybegin, &myend);
210 for ( i=mybegin; i<myend; ++i ) {
211#else
212 for ( i=0; i<m; ++i ) {
213#endif
214 iloc=rows[i];
215 for ( k=A->IA[iloc]; k<A->IA[iloc+1]; ++k ) {
216 j = A->JA[k];
217 if (mask[j]>=0) Aloc[i*n+mask[j]]=A->val[k];
218 } /* end for k */
219#ifdef _OPENMP
220 }
221 }
222#else
223 } /* enf for i */
224#endif
225
226#ifdef _OPENMP
227#pragma omp parallel for if(n>OPENMP_HOLDS) private(j)
228#endif
229 for ( j=0; j<n; ++j ) mask[cols[j]] = -1; // re-initialize mask
230
231 return FASP_SUCCESS;
232}
233
250static SHORT gentisquare_nomass (dCSRmat *A,
251 INT mm,
252 INT *Ii,
253 REAL *ima,
254 INT *mask)
255{
256 SHORT status = FASP_SUCCESS;
257
258 REAL *ms = (REAL *)fasp_mem_calloc(mm*mm,sizeof(REAL));
259
260 get_block(A,mm,mm,Ii,Ii,ms,mask);
261
262 status = invden(mm,ms,ima);
263
264 fasp_mem_free(ms); ms = NULL;
265
266 return status;
267}
268
288static SHORT getinonefull (INT **mat,
289 REAL **matval,
290 INT *lengths,
291 INT mm,
292 INT *Ii,
293 REAL *ima)
294{
295 INT tniz,i,j;
296
297#ifdef _OPENMP
298 // variables for OpenMP
299 INT myid, mybegin, myend;
300 INT nthreads = fasp_get_num_threads();
301#endif
302
303 tniz=lengths[1];
304
305#ifdef _OPENMP
306#pragma omp parallel for private(myid,mybegin,myend,i,j) if(mm>OPENMP_HOLDS)
307 for (myid=0; myid<nthreads; ++myid) {
308 fasp_get_start_end(myid, nthreads, mm, &mybegin, &myend);
309 for (i=mybegin; i<myend; ++i) {
310#else
311 for (i=0;i<mm;++i) {
312#endif
313 for (j=0;j<mm;++j) {
314 mat[0][tniz+i*mm+j]=Ii[i];
315 mat[1][tniz+i*mm+j]=Ii[j];
316 matval[0][tniz+i*mm+j]=ima[i*mm+j];
317 }
318#ifdef _OPENMP
319 }
320 }
321#else
322 }
323#endif
324 lengths[1]=tniz+mm*mm;
325
326 return FASP_SUCCESS;
327}
328
345static SHORT orderone (INT **mat,
346 REAL **matval,
347 INT *lengths)
348// lengths[0] for the number of rows
349// lengths[1] for the number of cols
350// lengths[2] for the number of nonzeros
351{
352 INT *rows[2],*cols[2],nns[2],tnizs[2];
353 REAL *vals[2];
354 SHORT status = FASP_SUCCESS;
355 INT tniz,i;
356
357 nns[0]=lengths[0];
358 nns[1]=lengths[1];
359 tnizs[0]=lengths[2];
360 tniz=lengths[2];
361
362 rows[0]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
363 cols[0]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
364 vals[0]=(REAL *)fasp_mem_calloc(tniz,sizeof(REAL));
365
366#ifdef _OPENMP
367#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
368#endif
369 for (i=0;i<tniz;++i) {
370 rows[0][i]=mat[0][i];
371 cols[0][i]=mat[1][i];
372 vals[0][i]=matval[0][i];
373 }
374
375 rows[1]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
376 cols[1]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
377 vals[1]=(REAL *)fasp_mem_calloc(tniz,sizeof(REAL));
378
379 fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
380
381 // all the nonzeros with same col are gathering together
382
383#ifdef _OPENMP
384#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
385#endif
386 for (i=0;i<tniz;++i) {
387 rows[0][i]=rows[1][i];
388 cols[0][i]=cols[1][i];
389 vals[0][i]=vals[1][i];
390 }
391 tnizs[1]=nns[0];
392 nns[0]=nns[1];
393 nns[1]=tnizs[1];
394 tnizs[1]=tnizs[0];
395 fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
396
397 // all the nonzeros with same col and row are gathering together
398#ifdef _OPENMP
399#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
400#endif
401 for (i=0;i<tniz;++i) {
402 rows[0][i]=rows[1][i];
403 cols[0][i]=cols[1][i];
404 vals[0][i]=vals[1][i];
405 }
406 tnizs[1]=nns[0];
407 nns[0]=nns[1];
408 nns[1]=tnizs[1];
409 tnizs[1]=tnizs[0];
410
411 tniz=tnizs[0];
412 for (i=0;i<tniz-1;++i) {
413 if (rows[0][i]==rows[0][i+1]&&cols[0][i]==cols[0][i+1]) {
414 vals[0][i+1]+=vals[0][i];
415 rows[0][i]=nns[0];
416 cols[0][i]=nns[1];
417 }
418 }
419 nns[0]=nns[0]+1;
420 nns[1]=nns[1]+1;
421
422 fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
423
424#ifdef _OPENMP
425#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
426#endif
427 for (i=0;i<tniz;++i) {
428 rows[0][i]=rows[1][i];
429 cols[0][i]=cols[1][i];
430 vals[0][i]=vals[1][i];
431 }
432 tnizs[1]=nns[0];
433 nns[0]=nns[1];
434 nns[1]=tnizs[1];
435 tnizs[1]=tnizs[0];
436
437 fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
438
439#ifdef _OPENMP
440#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
441#endif
442 for (i=0;i<tniz;++i) {
443 rows[0][i]=rows[1][i];
444 cols[0][i]=cols[1][i];
445 vals[0][i]=vals[1][i];
446 }
447 tnizs[1]=nns[0];
448 nns[0]=nns[1];
449 nns[1]=tnizs[1];
450 tnizs[1]=tnizs[0];
451
452 tniz=0;
453 for (i=0;i<tnizs[0];++i)
454 if (rows[0][i]<nns[0]-1) tniz++;
455
456#ifdef _OPENMP
457#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
458#endif
459 for (i=0;i<tniz;++i) {
460 mat[0][i]=rows[0][i];
461 mat[1][i]=cols[0][i];
462 matval[0][i]=vals[0][i];
463 }
464 nns[0]=nns[0]-1;
465 nns[1]=nns[1]-1;
466 lengths[0]=nns[0];
467 lengths[1]=nns[1];
468 lengths[2]=tniz;
469
470 fasp_mem_free(rows[0]); rows[0] = NULL;
471 fasp_mem_free(rows[1]); rows[1] = NULL;
472 fasp_mem_free(cols[0]); cols[0] = NULL;
473 fasp_mem_free(cols[1]); cols[1] = NULL;
474 fasp_mem_free(vals[0]); vals[0] = NULL;
475 fasp_mem_free(vals[1]); vals[1] = NULL;
476
477 return(status);
478}
479
511static SHORT genintval (dCSRmat *A,
512 INT **itmat,
513 REAL **itmatval,
514 INT ittniz,
515 INT *isol,
516 INT numiso,
517 INT nf,
518 INT nc)
519{
520 INT *Ii=NULL, *mask=NULL;
521 REAL *ima=NULL, *pex=NULL, **imas=NULL;
522 INT **mat=NULL;
523 REAL **matval;
524 INT lengths[3];
525 dCSRmat T;
526 INT tniz;
527 dvector sol, rhs;
528
529 INT mm,sum,i,j,k;
530 INT *iz,*izs,*izt,*izts;
531 SHORT status=FASP_SUCCESS;
532
533 mask=(INT *)fasp_mem_calloc(nf,sizeof(INT));
534 iz=(INT *)fasp_mem_calloc(nc,sizeof(INT));
535 izs=(INT *)fasp_mem_calloc(nc,sizeof(INT));
536 izt=(INT *)fasp_mem_calloc(nf,sizeof(INT));
537 izts=(INT *)fasp_mem_calloc(nf,sizeof(INT));
538
539 fasp_iarray_set(nf, mask, -1);
540
541 memset(iz, 0, sizeof(INT)*nc);
542
543#ifdef _OPENMP
544#pragma omp parallel for if(ittniz>OPENMP_HOLDS) private(i)
545#endif
546 for (i=0;i<ittniz;++i) iz[itmat[0][i]]++;
547
548 izs[0]=0;
549 for (i=1;i<nc;++i) izs[i]=izs[i-1]+iz[i-1];
550
551 sum = 0;
552#ifdef _OPENMP
553#pragma omp parallel for reduction(+:sum) if(nc>OPENMP_HOLDS) private(i)
554#endif
555 for (i=0;i<nc;++i) sum+=iz[i]*iz[i];
556
557 imas=(REAL **)fasp_mem_calloc(nc,sizeof(REAL *));
558
559 for (i=0;i<nc;++i) {
560 imas[i]=(REAL *)fasp_mem_calloc(iz[i]*iz[i],sizeof(REAL));
561 }
562
563 mat=(INT **)fasp_mem_calloc(2,sizeof(INT *));
564 mat[0]=(INT *)fasp_mem_calloc((sum+numiso),sizeof(INT));
565 mat[1]=(INT *)fasp_mem_calloc((sum+numiso),sizeof(INT));
566 matval=(REAL **)fasp_mem_calloc(1,sizeof(REAL *));
567 matval[0]=(REAL *)fasp_mem_calloc(sum+numiso,sizeof(REAL));
568
569 lengths[1]=0;
570
571 for (i=0;i<nc;++i) {
572
573 mm=iz[i];
574 Ii=(INT *)fasp_mem_realloc(Ii,mm*sizeof(INT));
575
576#ifdef _OPENMP
577#pragma omp parallel for if(mm>OPENMP_HOLDS) private(j)
578#endif
579 for (j=0;j<mm;++j) Ii[j]=itmat[1][izs[i]+j];
580
581 ima=(REAL *)fasp_mem_realloc(ima,mm*mm*sizeof(REAL));
582
583 gentisquare_nomass(A,mm,Ii,ima,mask);
584
585 getinonefull(mat,matval,lengths,mm,Ii,ima);
586
587#ifdef _OPENMP
588#pragma omp parallel for if(mm*mm>OPENMP_HOLDS) private(j)
589#endif
590 for (j=0;j<mm*mm;++j) imas[i][j]=ima[j];
591 }
592
593#ifdef _OPENMP
594#pragma omp parallel for if(numiso>OPENMP_HOLDS) private(i)
595#endif
596 for (i=0;i<numiso;++i) {
597 mat[0][sum+i]=isol[i];
598 mat[1][sum+i]=isol[i];
599 matval[0][sum+i]=1.0;
600 }
601
602 lengths[0]=nf;
603 lengths[2]=lengths[1]+numiso;
604 lengths[1]=nf;
605 orderone(mat,matval,lengths);
606 tniz=lengths[2];
607
608 sol.row=nf;
609 sol.val=(REAL*)fasp_mem_calloc(nf,sizeof(REAL));
610
611 memset(izt, 0, sizeof(INT)*nf);
612
613#ifdef _OPENMP
614#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(i)
615#endif
616 for (i=0;i<tniz;++i) izt[mat[0][i]]++;
617
618 T.IA=(INT*)fasp_mem_calloc((nf+1),sizeof(INT));
619
620 T.row=nf;
621 T.col=nf;
622 T.nnz=tniz;
623 T.IA[0]=0;
624 for (i=1;i<nf+1;++i) T.IA[i]=T.IA[i-1]+izt[i-1];
625
626 T.JA=(INT*)fasp_mem_calloc(tniz,sizeof(INT));
627
628#ifdef _OPENMP
629#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(j)
630#endif
631 for (j=0;j<tniz;++j) T.JA[j]=mat[1][j];
632
633 T.val=(REAL*)fasp_mem_calloc(tniz,sizeof(REAL));
634
635#ifdef _OPENMP
636#pragma omp parallel for if(tniz>OPENMP_HOLDS) private(j)
637#endif
638 for (j=0;j<tniz;++j) T.val[j]=matval[0][j];
639
640 rhs.val=(REAL*)fasp_mem_calloc(nf,sizeof(REAL));
641
642#ifdef _OPENMP
643#pragma omp parallel for if(nf>OPENMP_HOLDS) private(i)
644#endif
645 for (i=0;i<nf;++i) rhs.val[i]=1.0;
646 rhs.row=nf;
647
648 // setup preconditioner
649 dvector diag; fasp_dcsr_getdiag(0,&T,&diag);
650
651 precond pc;
652 pc.data = &diag;
653 pc.fct = fasp_precond_diag;
654
655 status = fasp_solver_dcsr_pcg(&T,&rhs,&sol,&pc,1e-3,1e-15,100,STOP_REL_RES,PRINT_NONE);
656
657 for (i=0;i<nc;++i) {
658 mm=iz[i];
659
660 ima=(REAL *)fasp_mem_realloc(ima,mm*mm*sizeof(REAL));
661
662 pex=(REAL *)fasp_mem_realloc(pex,mm*sizeof(REAL));
663
664 Ii=(INT *)fasp_mem_realloc(Ii,mm*sizeof(INT));
665
666#ifdef _OPENMP
667#pragma omp parallel for if(mm>OPENMP_HOLDS) private(j)
668#endif
669 for (j=0;j<mm;++j) Ii[j]=itmat[1][izs[i]+j];
670
671#ifdef _OPENMP
672#pragma omp parallel for if(mm*mm>OPENMP_HOLDS) private(j)
673#endif
674 for (j=0;j<mm*mm;++j) ima[j]=imas[i][j];
675
676#ifdef _OPENMP
677#pragma omp parallel for if(mm>OPENMP_HOLDS) private(k,j)
678#endif
679 for (k=0;k<mm;++k) {
680 for (pex[k]=j=0;j<mm;++j) pex[k]+=ima[k*mm+j]*sol.val[Ii[j]];
681 }
682#ifdef _OPENMP
683#pragma omp parallel for if(mm>OPENMP_HOLDS) private(j)
684#endif
685 for (j=0;j<mm;++j) itmatval[0][izs[i]+j]=pex[j];
686
687 }
688
689 fasp_mem_free(ima); ima = NULL;
690 fasp_mem_free(pex); pex = NULL;
691 fasp_mem_free(Ii); Ii = NULL;
692 fasp_mem_free(mask); mask = NULL;
693 fasp_mem_free(iz); iz = NULL;
694 fasp_mem_free(izs); izs = NULL;
695 fasp_mem_free(izt); izt = NULL;
696 fasp_mem_free(izts); izts = NULL;
697 fasp_mem_free(mat[0]); mat[0] = NULL;
698 fasp_mem_free(mat[1]); mat[1] = NULL;
699 fasp_mem_free(mat); mat = NULL;
700 fasp_mem_free(matval[0]); matval[0] = NULL;
701 fasp_mem_free(matval); matval = NULL;
702 for ( i=0; i<nc; ++i ) {fasp_mem_free(imas[i]); imas[i] = NULL;}
703 fasp_mem_free(imas); imas = NULL;
704
705 fasp_dcsr_free(&T);
706 fasp_dvec_free(&rhs);
707 fasp_dvec_free(&sol);
708 fasp_dvec_free(&diag);
709
710 return status;
711}
712
727static SHORT getiteval (dCSRmat *A,
728 dCSRmat *it)
729{
730 INT nf,nc,ittniz;
731 INT *itmat[2];
732 REAL **itmatval;
733 INT *rows[2],*cols[2];
734 REAL *vals[2];
735 INT nns[2],tnizs[2];
736 INT i,j,numiso;
737 INT *isol;
738 SHORT status = FASP_SUCCESS;
739
740 nf=A->row;
741 nc=it->col;
742 ittniz=it->IA[nf];
743
744 itmat[0]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
745 itmat[1]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
746 itmatval=(REAL **)fasp_mem_calloc(1,sizeof(REAL *));
747 itmatval[0]=(REAL *)fasp_mem_calloc(ittniz,sizeof(REAL));
748 isol=(INT *)fasp_mem_calloc(nf,sizeof(INT));
749
750 numiso=0;
751 for (i=0;i<nf;++i) {
752 if (it->IA[i]==it->IA[i+1]) {
753 isol[numiso]=i;
754 numiso++;
755 }
756 }
757
758#ifdef _OPENMP
759#pragma omp parallel for if(nf>OPENMP_HOLDS) private(i,j)
760#endif
761 for (i=0;i<nf;++i) {
762 for (j=it->IA[i];j<it->IA[i+1];++j) itmat[0][j]=i;
763 }
764
765#ifdef _OPENMP
766#pragma omp parallel for if(ittniz>OPENMP_HOLDS) private(j)
767#endif
768 for (j=0;j<ittniz;++j) {
769 itmat[1][j]=it->JA[j];
770 itmatval[0][j]=it->val[j];
771 }
772
773 rows[0]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
774 cols[0]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
775 vals[0]=(REAL *)fasp_mem_calloc(ittniz,sizeof(REAL));
776
777#ifdef _OPENMP
778#pragma omp parallel for if(ittniz>OPENMP_HOLDS) private(i)
779#endif
780 for (i=0;i<ittniz;++i) {
781 rows[0][i]=itmat[0][i];
782 cols[0][i]=itmat[1][i];
783 vals[0][i]=itmat[0][i];
784 }
785
786 nns[0]=nf;
787 nns[1]=nc;
788 tnizs[0]=ittniz;
789
790 rows[1]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
791 cols[1]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
792 vals[1]=(REAL *)fasp_mem_calloc(ittniz,sizeof(REAL));
793
794 fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
795
796#ifdef _OPENMP
797#pragma omp parallel for if(ittniz>OPENMP_HOLDS) private(i)
798#endif
799 for (i=0;i<ittniz;++i) {
800 itmat[0][i]=rows[1][i];
801 itmat[1][i]=cols[1][i];
802 itmatval[0][i]=vals[1][i];
803 }
804 genintval(A,itmat,itmatval,ittniz,isol,numiso,nf,nc);
805
806#ifdef _OPENMP
807#pragma omp parallel for if(ittniz>OPENMP_HOLDS) private(i)
808#endif
809 for (i=0;i<ittniz;++i) {
810 rows[0][i]=itmat[0][i];
811 cols[0][i]=itmat[1][i];
812 vals[0][i]=itmatval[0][i];
813 }
814 nns[0]=nc;
815 nns[1]=nf;
816 tnizs[0]=ittniz;
817
818 fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
819
820#ifdef _OPENMP
821#pragma omp parallel for if(ittniz>OPENMP_HOLDS) private(i)
822#endif
823 for (i=0;i<ittniz;++i) it->val[i]=vals[1][i];
824
825 fasp_mem_free(isol); isol = NULL;
826 fasp_mem_free(itmat[0]); itmat[0] = NULL;
827 fasp_mem_free(itmat[1]); itmat[1] = NULL;
828 fasp_mem_free(itmatval[0]); itmatval[0] = NULL;
829 fasp_mem_free(itmatval); itmatval = NULL;
830 fasp_mem_free(rows[0]); rows[0] = NULL;
831 fasp_mem_free(rows[1]); rows[1] = NULL;
832 fasp_mem_free(cols[0]); cols[0] = NULL;
833 fasp_mem_free(cols[1]); cols[1] = NULL;
834 fasp_mem_free(vals[0]); vals[0] = NULL;
835 fasp_mem_free(vals[1]); vals[1] = NULL;
836
837 return status;
838}
839
840/*---------------------------------*/
841/*-- End of File --*/
842/*---------------------------------*/
fasp_iarray_set
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:98
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_realloc
void * fasp_mem_realloc(void *oldmem, const LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: AuxMemory.c:113
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_get_start_end
void fasp_get_start_end(const INT procid, const INT nprocs, const INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: AuxThreads.c:93
fasp_dvec_free
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
Definition: AuxVector.c:145
fasp_smat_lu_decomp
SHORT fasp_smat_lu_decomp(REAL *A, INT pivot[], const INT n)
LU decomposition of A using Doolittle's method.
Definition: BlaSmallMatLU.c:52
fasp_smat_lu_solve
SHORT fasp_smat_lu_solve(const REAL *A, REAL b[], const INT pivot[], REAL x[], const INT n)
Solving Ax=b using LU decomposition.
Definition: BlaSmallMatLU.c:124
fasp_dcsr_free
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:184
fasp_dcsr_transpose
void fasp_dcsr_transpose(INT *row[2], INT *col[2], REAL *val[2], INT *nn, INT *tniz)
Transpose of a dCSRmat matrix.
Definition: BlaSparseCSR.c:1037
fasp_dcsr_getdiag
void fasp_dcsr_getdiag(INT n, const dCSRmat *A, dvector *diag)
Get first n diagonal entries of a CSR matrix A.
Definition: BlaSparseCSR.c:537
fasp_solver_dcsr_pcg
INT fasp_solver_dcsr_pcg(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: KryPcg.c:96
fasp_amg_interp_em
void fasp_amg_interp_em(dCSRmat *A, ivector *vertices, dCSRmat *P, AMG_param *param)
Energy-min interpolation.
Definition: PreAMGInterpEM.c:63
fasp_precond_diag
void fasp_precond_diag(REAL *r, REAL *z, void *data)
Diagonal preconditioner z=inv(D)*r.
Definition: PreCSR.c:172
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
INT
#define INT
Definition: fasp.h:72
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
STOP_REL_RES
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:132
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
AMG_param
Parameters for AMG methods.
Definition: fasp.h:455
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
dvector::row
INT row
number of rows
Definition: fasp.h:357
ivector
Vector with n entries of INT type.
Definition: fasp.h:368
ivector::row
INT row
number of rows
Definition: fasp.h:371
ivector::val
INT * val
actual vector entries
Definition: fasp.h:374
precond
Preconditioner data and action.
Definition: fasp.h:1095
precond::data
void * data
data for preconditioner, void pointer
Definition: fasp.h:1098
precond::fct
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1101