Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
interpolation_em.c
Go to the documentation of this file.
1 
10 #include <math.h>
11 #include <time.h>
12 
13 #ifdef _OPENMP
14 #include <omp.h>
15 #endif
16 
17 #include "fasp.h"
18 #include "fasp_functs.h"
19 
20 static SHORT getiteval(dCSRmat *, dCSRmat *);
21 static SHORT invden(INT, REAL *, REAL *);
22 static SHORT get_block (dCSRmat *, INT, INT, INT *, INT *, REAL *, INT *);
23 static SHORT gentisquare_nomass(dCSRmat *, INT, INT *, REAL *, INT *);
24 static SHORT getinonefull(INT **, REAL **, INT *, INT, INT *, REAL *);
25 static SHORT orderone(INT **, REAL **, INT *);
26 static SHORT genintval(dCSRmat *, INT **, REAL **, INT, INT *, INT, INT, INT);
27 
28 /*---------------------------------*/
29 /*-- Public Functions --*/
30 /*---------------------------------*/
31 
50  ivector *vertices,
51  dCSRmat *P,
52  AMG_param *param)
53 {
54  INT *vec = vertices->val;
55  INT *CoarseIndex = (INT *)fasp_mem_calloc(vertices->row, sizeof(INT));
56  INT i, j, index;
57 
58  // generate indices for C-points
59  for ( index = i = 0; i < vertices->row; ++i ) {
60  if ( vec[i] == 1 ) {
61  CoarseIndex[i] = index;
62  index++;
63  }
64  }
65 
66 #ifdef _OPENMP
67 #pragma omp parallel for private(i,j) if(P->nnz>OPENMP_HOLDS)
68 #endif
69  for ( i = 0; i < P->nnz; ++i ) {
70  j = P->JA[i];
71  P->JA[i] = CoarseIndex[j];
72  }
73 
74  // clean up memory
75  fasp_mem_free(CoarseIndex);
76 
77  // main part
78  getiteval(A, P);
79 }
80 
81 /*---------------------------------*/
82 /*-- Private Functions --*/
83 /*---------------------------------*/
84 
103 static SHORT invden (INT nn,
104  REAL *mat,
105  REAL *invmat)
106 {
107  INT i,j;
108  SHORT status = FASP_SUCCESS;
109 
110 #ifdef _OPENMP
111  // variables for OpenMP
112  INT myid, mybegin, myend;
113  INT nthreads = FASP_GET_NUM_THREADS();
114 #endif
115 
116  INT *pivot=(INT *)fasp_mem_calloc(nn,sizeof(INT));
117  REAL *rhs=(REAL *)fasp_mem_calloc(nn,sizeof(REAL));
118  REAL *sol=(REAL *)fasp_mem_calloc(nn,sizeof(REAL));
119 
120  fasp_smat_lu_decomp(mat,pivot,nn);
121 
122 #ifdef _OPENMP
123 #pragma omp parallel for private(myid,mybegin,myend,i,j) if(nn>OPENMP_HOLDS)
124  for (myid=0; myid<nthreads; ++myid) {
125  FASP_GET_START_END(myid, nthreads, nn, &mybegin, &myend);
126  for (i=mybegin; i<myend; ++i) {
127 #else
128  for (i=0;i<nn;++i) {
129 #endif
130  for (j=0;j<nn;++j) rhs[j]=0.;
131  rhs[i]=1.;
132  fasp_smat_lu_solve(mat,rhs,pivot,sol,nn);
133  for (j=0;j<nn;++j) invmat[i*nn+j]=sol[j];
134 #ifdef _OPENMP
135  }
136  }
137 #else
138  }
139 #endif
140 
141  fasp_mem_free(pivot);
142  fasp_mem_free(rhs);
143  fasp_mem_free(sol);
144 
145  return status;
146 }
147 
169 static SHORT get_block (dCSRmat *A,
170  INT m,
171  INT n,
172  INT *rows,
173  INT *cols,
174  REAL *Aloc,
175  INT *mask)
176 {
177  INT i, j, k, iloc;
178 
179 #ifdef _OPENMP
180  // variables for OpenMP
181  INT myid, mybegin, myend;
182  INT nthreads = FASP_GET_NUM_THREADS();
183 #endif
184 
185 #if 0
186  for ( i=0; i<m; ++i ) {
187  for ( j=0; j<n; ++j ) {
188  Aloc[i*n+j] = 0.0; // initialize Aloc
189  }
190  }
191 #endif
192  memset(Aloc, 0x0, sizeof(REAL)*m*n);
193 
194 #ifdef _OPENMP
195 #pragma omp parallel for if(n>OPENMP_HOLDS)
196 #endif
197  for ( j=0; j<n; ++j ) {
198  mask[cols[j]] = j; // initialize mask, mask stores C indices 0,1,...
199  }
200 
201 #ifdef _OPENMP
202 #pragma omp parallel for private(myid,mybegin,myend,i,j,k,iloc) if(m>OPENMP_HOLDS)
203  for ( myid=0; myid<nthreads; ++myid ) {
204  FASP_GET_START_END(myid, nthreads, m, &mybegin, &myend);
205  for ( i=mybegin; i<myend; ++i ) {
206 #else
207  for ( i=0; i<m; ++i ) {
208 #endif
209  iloc=rows[i];
210  for ( k=A->IA[iloc]; k<A->IA[iloc+1]; ++k ) {
211  j = A->JA[k];
212  if (mask[j]>=0) Aloc[i*n+mask[j]]=A->val[k];
213  } /* end for k */
214 #ifdef _OPENMP
215  }
216  }
217 #else
218  } /* enf for i */
219 #endif
220 
221 #ifdef _OPENMP
222 #pragma omp parallel for if(n>OPENMP_HOLDS)
223 #endif
224  for ( j=0; j<n; ++j ) mask[cols[j]] = -1; // re-initialize mask
225 
226  return FASP_SUCCESS;
227 }
228 
245 static SHORT gentisquare_nomass (dCSRmat *A,
246  INT mm,
247  INT *Ii,
248  REAL *ima,
249  INT *mask)
250 {
251  SHORT status = FASP_SUCCESS;
252 
253  REAL *ms = (REAL *)fasp_mem_calloc(mm*mm,sizeof(REAL));
254 
255  status = get_block(A,mm,mm,Ii,Ii,ms,mask);
256  status = invden(mm,ms,ima);
257 
258  fasp_mem_free(ms);
259  return status;
260 }
261 
281 static SHORT getinonefull (INT **mat,
282  REAL **matval,
283  INT *lengths,
284  INT mm,
285  INT *Ii,
286  REAL *ima)
287 {
288  INT tniz,i,j;
289 
290 #ifdef _OPENMP
291  // variables for OpenMP
292  INT myid, mybegin, myend;
293  INT nthreads = FASP_GET_NUM_THREADS();
294 #endif
295 
296  tniz=lengths[1];
297 
298 #ifdef _OPENMP
299 #pragma omp parallel for private(myid,mybegin,myend,i,j) if(mm>OPENMP_HOLDS)
300  for (myid=0; myid<nthreads; ++myid) {
301  FASP_GET_START_END(myid, nthreads, mm, &mybegin, &myend);
302  for (i=mybegin; i<myend; ++i) {
303 #else
304  for (i=0;i<mm;++i) {
305 #endif
306  for (j=0;j<mm;++j) {
307  mat[0][tniz+i*mm+j]=Ii[i];
308  mat[1][tniz+i*mm+j]=Ii[j];
309  matval[0][tniz+i*mm+j]=ima[i*mm+j];
310  }
311 #ifdef _OPENMP
312  }
313  }
314 #else
315  }
316 #endif
317  lengths[1]=tniz+mm*mm;
318 
319  return FASP_SUCCESS;
320 }
321 
338 static SHORT orderone (INT **mat,
339  REAL **matval,
340  INT *lengths)
341 // lengths[0] for the number of rows
342 // lengths[1] for the number of cols
343 // lengths[2] for the number of nonzeros
344 {
345  INT *rows[2],*cols[2],nns[2],tnizs[2];
346  REAL *vals[2];
347  SHORT status = FASP_SUCCESS;
348  INT tniz,i;
349 
350  nns[0]=lengths[0];
351  nns[1]=lengths[1];
352  tnizs[0]=lengths[2];
353  tniz=lengths[2];
354 
355  rows[0]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
356  cols[0]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
357  vals[0]=(REAL *)fasp_mem_calloc(tniz,sizeof(REAL));
358 
359 #ifdef _OPENMP
360 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
361 #endif
362  for (i=0;i<tniz;++i) {
363  rows[0][i]=mat[0][i];
364  cols[0][i]=mat[1][i];
365  vals[0][i]=matval[0][i];
366  }
367 
368  rows[1]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
369  cols[1]=(INT *)fasp_mem_calloc(tniz,sizeof(INT));
370  vals[1]=(REAL *)fasp_mem_calloc(tniz,sizeof(REAL));
371 
372  fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
373 
374  // all the nonzeros with same col are gathering together
375 
376 #ifdef _OPENMP
377 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
378 #endif
379  for (i=0;i<tniz;++i) {
380  rows[0][i]=rows[1][i];
381  cols[0][i]=cols[1][i];
382  vals[0][i]=vals[1][i];
383  }
384  tnizs[1]=nns[0];
385  nns[0]=nns[1];
386  nns[1]=tnizs[1];
387  tnizs[1]=tnizs[0];
388  fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
389 
390  // all the nonzeros with same col and row are gathering together
391 
392 #ifdef _OPENMP
393 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
394 #endif
395  for (i=0;i<tniz;++i) {
396  rows[0][i]=rows[1][i];
397  cols[0][i]=cols[1][i];
398  vals[0][i]=vals[1][i];
399  }
400  tnizs[1]=nns[0];
401  nns[0]=nns[1];
402  nns[1]=tnizs[1];
403  tnizs[1]=tnizs[0];
404 
405  tniz=tnizs[0];
406  for (i=0;i<tniz-1;++i) {
407  if (rows[0][i]==rows[0][i+1]&&cols[0][i]==cols[0][i+1]) {
408  vals[0][i+1]+=vals[0][i];
409  rows[0][i]=nns[0];
410  cols[0][i]=nns[1];
411  }
412  }
413  nns[0]=nns[0]+1;
414  nns[1]=nns[1]+1;
415 
416  fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
417 
418 #ifdef _OPENMP
419 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
420 #endif
421  for (i=0;i<tniz;++i) {
422  rows[0][i]=rows[1][i];
423  cols[0][i]=cols[1][i];
424  vals[0][i]=vals[1][i];
425  }
426  tnizs[1]=nns[0];
427  nns[0]=nns[1];
428  nns[1]=tnizs[1];
429  tnizs[1]=tnizs[0];
430 
431  fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
432 
433 #ifdef _OPENMP
434 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
435 #endif
436  for (i=0;i<tniz;++i) {
437  rows[0][i]=rows[1][i];
438  cols[0][i]=cols[1][i];
439  vals[0][i]=vals[1][i];
440  }
441  tnizs[1]=nns[0];
442  nns[0]=nns[1];
443  nns[1]=tnizs[1];
444  tnizs[1]=tnizs[0];
445 
446  tniz=0;
447  for (i=0;i<tnizs[0];++i)
448  if (rows[0][i]<nns[0]-1) tniz++;
449 
450 #ifdef _OPENMP
451 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
452 #endif
453  for (i=0;i<tniz;++i) {
454  mat[0][i]=rows[0][i];
455  mat[1][i]=cols[0][i];
456  matval[0][i]=vals[0][i];
457  }
458  nns[0]=nns[0]-1;
459  nns[1]=nns[1]-1;
460  lengths[0]=nns[0];
461  lengths[1]=nns[1];
462  lengths[2]=tniz;
463 
464  fasp_mem_free(rows[0]);
465  fasp_mem_free(rows[1]);
466  fasp_mem_free(cols[0]);
467  fasp_mem_free(cols[1]);
468  fasp_mem_free(vals[0]);
469  fasp_mem_free(vals[1]);
470 
471  return(status);
472 }
473 
474 
506 static SHORT genintval (dCSRmat *A,
507  INT **itmat,
508  REAL **itmatval,
509  INT ittniz,
510  INT *isol,
511  INT numiso,
512  INT nf,
513  INT nc)
514 {
515  INT *Ii=NULL, *mask=NULL;
516  REAL *ima=NULL, *pex=NULL, **imas=NULL;
517  INT **mat=NULL;
518  REAL **matval;
519  INT lengths[3];
520  dCSRmat T;
521  INT tniz;
522  dvector sol, rhs;
523 
524  INT mm,sum,i,j,k;
525  INT *iz,*izs,*izt,*izts;
526  SHORT status=FASP_SUCCESS;
527 
528  itsolver_param itparam;
529  itparam.print_level = PRINT_NONE;
530  itparam.itsolver_type = SOLVER_CG;
531  itparam.stop_type = STOP_REL_RES;
532  itparam.tol = 1e-3;
533  itparam.maxit = 100;
534  itparam.restart = 100;
535 
536  mask=(INT *)fasp_mem_calloc(nf,sizeof(INT));
537  iz=(INT *)fasp_mem_calloc(nc,sizeof(INT));
538  izs=(INT *)fasp_mem_calloc(nc,sizeof(INT));
539  izt=(INT *)fasp_mem_calloc(nf,sizeof(INT));
540  izts=(INT *)fasp_mem_calloc(nf,sizeof(INT));
541 
542  fasp_iarray_set(nf, mask, -1);
543 
544  memset(iz, 0, sizeof(INT)*nc);
545 
546 #ifdef _OPENMP
547 #pragma omp parallel for if(ittniz>OPENMP_HOLDS)
548 #endif
549  for (i=0;i<ittniz;++i) iz[itmat[0][i]]++;
550 
551  izs[0]=0;
552  for (i=1;i<nc;++i) izs[i]=izs[i-1]+iz[i-1];
553 
554  sum = 0;
555 #ifdef _OPENMP
556 #pragma omp parallel for reduction(+:sum) if(nc>OPENMP_HOLDS)
557 #endif
558  for (i=0;i<nc;++i) sum+=iz[i]*iz[i];
559 
560  imas=(REAL **)fasp_mem_calloc(nc,sizeof(REAL *));
561 
562  for (i=0;i<nc;++i) {
563  imas[i]=(REAL *)fasp_mem_calloc(iz[i]*iz[i],sizeof(REAL));
564  }
565 
566  mat=(INT **)fasp_mem_calloc(2,sizeof(INT *));
567  mat[0]=(INT *)fasp_mem_calloc((sum+numiso),sizeof(INT));
568  mat[1]=(INT *)fasp_mem_calloc((sum+numiso),sizeof(INT));
569  matval=(REAL **)fasp_mem_calloc(1,sizeof(REAL *));
570  matval[0]=(REAL *)fasp_mem_calloc(sum+numiso,sizeof(REAL));
571 
572  lengths[1]=0;
573 
574  for (i=0;i<nc;++i) {
575 
576  mm=iz[i];
577  Ii=(INT *)fasp_mem_realloc(Ii,mm*sizeof(INT));
578 
579 #ifdef _OPENMP
580 #pragma omp parallel for if(mm>OPENMP_HOLDS)
581 #endif
582  for (j=0;j<mm;++j) Ii[j]=itmat[1][izs[i]+j];
583 
584  ima=(REAL *)fasp_mem_realloc(ima,mm*mm*sizeof(REAL));
585 
586  gentisquare_nomass(A,mm,Ii,ima,mask);
587 
588  getinonefull(mat,matval,lengths,mm,Ii,ima);
589 
590 #ifdef _OPENMP
591 #pragma omp parallel for if(mm*mm>OPENMP_HOLDS)
592 #endif
593  for (j=0;j<mm*mm;++j) imas[i][j]=ima[j];
594  }
595 
596 #ifdef _OPENMP
597 #pragma omp parallel for if(numiso>OPENMP_HOLDS) private(i)
598 #endif
599  for (i=0;i<numiso;++i) {
600  mat[0][sum+i]=isol[i];
601  mat[1][sum+i]=isol[i];
602  matval[0][sum+i]=1.0;
603  }
604 
605  lengths[0]=nf;
606  lengths[2]=lengths[1]+numiso;
607  lengths[1]=nf;
608  orderone(mat,matval,lengths);
609  tniz=lengths[2];
610 
611  sol.row=nf;
612  sol.val=(REAL*)fasp_mem_calloc(nf,sizeof(REAL));
613 
614  memset(izt, 0, sizeof(INT)*nf);
615 
616 #ifdef _OPENMP
617 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
618 #endif
619  for (i=0;i<tniz;++i) izt[mat[0][i]]++;
620 
621  T.IA=(INT*)fasp_mem_calloc((nf+1),sizeof(INT));
622 
623  T.row=nf;
624  T.col=nf;
625  T.nnz=tniz;
626  T.IA[0]=0;
627  for (i=1;i<nf+1;++i) T.IA[i]=T.IA[i-1]+izt[i-1];
628 
629  T.JA=(INT*)fasp_mem_calloc(tniz,sizeof(INT));
630 
631 #ifdef _OPENMP
632 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
633 #endif
634  for (j=0;j<tniz;++j) T.JA[j]=mat[1][j];
635 
636  T.val=(REAL*)fasp_mem_calloc(tniz,sizeof(REAL));
637 
638 #ifdef _OPENMP
639 #pragma omp parallel for if(tniz>OPENMP_HOLDS)
640 #endif
641  for (j=0;j<tniz;++j) T.val[j]=matval[0][j];
642 
643  rhs.val=(REAL*)fasp_mem_calloc(nf,sizeof(REAL));
644 
645 #ifdef _OPENMP
646 #pragma omp parallel for if(nf>OPENMP_HOLDS)
647 #endif
648  for (i=0;i<nf;++i) rhs.val[i]=1.0;
649  rhs.row=nf;
650 
651  fasp_solver_dcsr_krylov_diag(&T,&rhs,&sol,&itparam);
652 
653  for (i=0;i<nc;++i) {
654  mm=iz[i];
655 
656  ima=(REAL *)fasp_mem_realloc(ima,mm*mm*sizeof(REAL));
657 
658  pex=(REAL *)fasp_mem_realloc(pex,mm*sizeof(REAL));
659 
660  Ii=(INT *)fasp_mem_realloc(Ii,mm*sizeof(INT));
661 
662 #ifdef _OPENMP
663 #pragma omp parallel for if(mm>OPENMP_HOLDS)
664 #endif
665  for (j=0;j<mm;++j) Ii[j]=itmat[1][izs[i]+j];
666 
667 #ifdef _OPENMP
668 #pragma omp parallel for if(mm*mm>OPENMP_HOLDS)
669 #endif
670  for (j=0;j<mm*mm;++j) ima[j]=imas[i][j];
671 
672 #ifdef _OPENMP
673 #pragma omp parallel for if(mm>OPENMP_HOLDS) private(k,j)
674 #endif
675  for (k=0;k<mm;++k) {
676  for (pex[k]=j=0;j<mm;++j) pex[k]+=ima[k*mm+j]*sol.val[Ii[j]];
677  }
678 #ifdef _OPENMP
679 #pragma omp parallel for if(mm>OPENMP_HOLDS)
680 #endif
681  for (j=0;j<mm;++j) itmatval[0][izs[i]+j]=pex[j];
682 
683  }
684 
685  fasp_mem_free(ima);
686  fasp_mem_free(pex);
687  fasp_mem_free(Ii);
688  fasp_mem_free(mask);
689  fasp_mem_free(iz);
690  fasp_mem_free(izs);
691  fasp_mem_free(izt);
692  fasp_mem_free(izts);
693  fasp_mem_free(mat[0]);
694  fasp_mem_free(mat[1]);
695  fasp_mem_free(mat);
696  fasp_mem_free(matval[0]);
697  fasp_mem_free(matval);
698  for (i=0;i<nc;++i) fasp_mem_free(imas[i]);
699  fasp_mem_free(imas);
700  fasp_dcsr_free(&T);
701  fasp_dvec_free(&rhs);
702  fasp_dvec_free(&sol);
703 
704  return status;
705 }
706 
721 static SHORT getiteval (dCSRmat *A,
722  dCSRmat *it)
723 {
724  INT nf,nc,ittniz;
725  INT *itmat[2];
726  REAL **itmatval;
727  INT *rows[2],*cols[2];
728  REAL *vals[2];
729  INT nns[2],tnizs[2];
730  INT i,j,numiso;
731  INT *isol;
732  SHORT status = FASP_SUCCESS;
733 
734  nf=A->row;
735  nc=it->col;
736  ittniz=it->IA[nf];
737 
738  itmat[0]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
739  itmat[1]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
740  itmatval=(REAL **)fasp_mem_calloc(1,sizeof(REAL *));
741  itmatval[0]=(REAL *)fasp_mem_calloc(ittniz,sizeof(REAL));
742  isol=(INT *)fasp_mem_calloc(nf,sizeof(INT));
743 
744  numiso=0;
745  for (i=0;i<nf;++i) {
746  if (it->IA[i]==it->IA[i+1]) {
747  isol[numiso]=i;
748  numiso++;
749  }
750  }
751 
752 #ifdef _OPENMP
753 #pragma omp parallel for if(nf>OPENMP_HOLDS) private(i,j)
754 #endif
755  for (i=0;i<nf;++i) {
756  for (j=it->IA[i];j<it->IA[i+1];++j) itmat[0][j]=i;
757  }
758 #if 0
759  for (j=0;j<ittniz;++j) itmat[1][j]=it->JA[j];
760  for (j=0;j<ittniz;++j) itmatval[0][j]=it->val[j];
761 #endif
762 
763 #ifdef _OPENMP
764 #pragma omp parallel for if(ittniz>OPENMP_HOLDS)
765 #endif
766  for (j=0;j<ittniz;++j) {
767  itmat[1][j]=it->JA[j];
768  itmatval[0][j]=it->val[j];
769  }
770 
771  rows[0]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
772  cols[0]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
773  vals[0]=(REAL *)fasp_mem_calloc(ittniz,sizeof(REAL));
774 
775 #ifdef _OPENMP
776 #pragma omp parallel for if(ittniz>OPENMP_HOLDS)
777 #endif
778  for (i=0;i<ittniz;++i) {
779  rows[0][i]=itmat[0][i];
780  cols[0][i]=itmat[1][i];
781  vals[0][i]=itmat[0][i];
782  }
783 
784  nns[0]=nf;
785  nns[1]=nc;
786  tnizs[0]=ittniz;
787 
788  rows[1]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
789  cols[1]=(INT *)fasp_mem_calloc(ittniz,sizeof(INT));
790  vals[1]=(REAL *)fasp_mem_calloc(ittniz,sizeof(REAL));
791 
792  fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
793 
794 #ifdef _OPENMP
795 #pragma omp parallel for if(ittniz>OPENMP_HOLDS)
796 #endif
797  for (i=0;i<ittniz;++i) {
798  itmat[0][i]=rows[1][i];
799  itmat[1][i]=cols[1][i];
800  itmatval[0][i]=vals[1][i];
801  }
802  genintval(A,itmat,itmatval,ittniz,isol,numiso,nf,nc);
803 
804 #ifdef _OPENMP
805 #pragma omp parallel for if(ittniz>OPENMP_HOLDS)
806 #endif
807  for (i=0;i<ittniz;++i) {
808  rows[0][i]=itmat[0][i];
809  cols[0][i]=itmat[1][i];
810  vals[0][i]=itmatval[0][i];
811  }
812  nns[0]=nc;
813  nns[1]=nf;
814  tnizs[0]=ittniz;
815 
816  fasp_dcsr_transpose(rows,cols,vals,nns,tnizs);
817 
818 #ifdef _OPENMP
819 #pragma omp parallel for if(ittniz>OPENMP_HOLDS)
820 #endif
821  for (i=0;i<ittniz;++i) it->val[i]=vals[1][i];
822 
823  fasp_mem_free(isol);
824  fasp_mem_free(itmat[0]);
825  fasp_mem_free(itmat[1]);
826  fasp_mem_free(itmatval[0]);
827  fasp_mem_free(itmatval);
828  fasp_mem_free(rows[0]);
829  fasp_mem_free(rows[1]);
830  fasp_mem_free(cols[0]);
831  fasp_mem_free(cols[1]);
832  fasp_mem_free(vals[0]);
833  fasp_mem_free(vals[1]);
834 
835  return status;
836 }
837 
838 /*---------------------------------*/
839 /*-- End of File --*/
840 /*---------------------------------*/
Parameters for AMG solver.
Definition: fasp.h:583
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
void fasp_amg_interp_em(dCSRmat *A, ivector *vertices, dCSRmat *P, AMG_param *param)
Energy-min interpolation.
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
Definition: vec.c:139
REAL * val
nonzero entries of A
Definition: fasp.h:166
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
Vector with n entries of REAL type.
Definition: fasp.h:342
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
INT restart
Definition: fasp.h:1112
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
SHORT stop_type
Definition: fasp.h:1109
#define OPENMP_HOLDS
Definition: fasp_const.h:248
INT nnz
number of nonzero entries
Definition: fasp.h:157
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
void * fasp_mem_realloc(void *oldmem, LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: memory.c:110
#define SOLVER_CG
Definition: fasp_const.h:103
SHORT fasp_smat_lu_decomp(REAL *A, INT pivot[], const INT n)
LU decomposition of A usind Doolittle's method.
Definition: lu.c:46
INT * val
actual vector entries
Definition: fasp.h:362
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
INT row
row number of matrix A, m
Definition: fasp.h:151
Parameters passed to iterative solvers.
Definition: fasp.h:1105
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:131
INT row
number of rows
Definition: fasp.h:345
REAL tol
Definition: fasp.h:1111
SHORT print_level
Definition: fasp.h:1113
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT fasp_solver_dcsr_krylov_diag(dCSRmat *A, dvector *b, dvector *x, itsolver_param *itparam)
Solve Ax=b by diagonal preconditioned Krylov methods.
Definition: itsolver_csr.c:193
INT row
number of rows
Definition: fasp.h:359
SHORT fasp_smat_lu_solve(REAL *A, REAL b[], INT pivot[], REAL x[], const INT n)
Solving Ax=b using LU decomposition.
Definition: lu.c:117
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: array.c:107
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:166
SHORT itsolver_type
Definition: fasp.h:1107