Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : #include <petsc/private/vecimpl.h>
      12             : 
      13             : #if defined(__WITH_MPI__)
      14             : #define __SUF__(A) A##_MPI
      15             : #else
      16             : #define __SUF__(A) A##_Seq
      17             : #endif
      18             : #define __QUOTEME_(x) #x
      19             : #define __QUOTEME(x) __QUOTEME_(x)
      20             : #define __SUF_C__(A) __QUOTEME(__SUF__(A))
      21             : 
      22       74291 : static PetscErrorCode __SUF__(VecDot_Comp)(Vec a,Vec b,PetscScalar *z)
      23             : {
      24       74291 :   PetscScalar    sum = 0.0,work;
      25       74291 :   PetscInt       i;
      26       74291 :   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
      27             : 
      28       74291 :   PetscFunctionBegin;
      29       74291 :   SlepcValidVecComp(a,1);
      30       74291 :   SlepcValidVecComp(b,2);
      31       74291 :   if (as->x[0]->ops->dot_local) {
      32      148585 :     for (i=0,sum=0.0;i<as->n->n;i++) {
      33       74294 :       PetscUseTypeMethod(as->x[i],dot_local,bs->x[i],&work);
      34       74294 :       sum += work;
      35             :     }
      36             : #if defined(__WITH_MPI__)
      37       74291 :     work = sum;
      38      222873 :     PetscCallMPI(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
      39             : #endif
      40             :   } else {
      41           0 :     for (i=0,sum=0.0;i<as->n->n;i++) {
      42           0 :       PetscCall(VecDot(as->x[i],bs->x[i],&work));
      43           0 :       sum += work;
      44             :     }
      45             :   }
      46       74291 :   *z = sum;
      47       74291 :   PetscFunctionReturn(PETSC_SUCCESS);
      48             : }
      49             : 
      50       54285 : static PetscErrorCode __SUF__(VecMDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
      51             : {
      52       54285 :   PetscScalar    *work,*work0,*r;
      53       54285 :   Vec_Comp       *as = (Vec_Comp*)a->data;
      54       54285 :   Vec            *bx;
      55       54285 :   PetscInt       i,j;
      56             : 
      57       54285 :   PetscFunctionBegin;
      58       54285 :   SlepcValidVecComp(a,1);
      59      185690 :   for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
      60             : 
      61       54285 :   if (as->n->n == 0) {
      62           0 :     *z = 0;
      63           0 :     PetscFunctionReturn(PETSC_SUCCESS);
      64             :   }
      65             : 
      66       54285 :   PetscCall(PetscMalloc2(n,&work0,n,&bx));
      67             : 
      68             : #if defined(__WITH_MPI__)
      69       54285 :   if (as->x[0]->ops->mdot_local) {
      70       54285 :     r = work0; work = z;
      71             :   } else
      72             : #endif
      73             :   {
      74           0 :     r = z; work = work0;
      75             :   }
      76             : 
      77             :   /* z[i] <- a.x' * b[i].x */
      78      185690 :   for (i=0;i<n;i++) r[i] = 0.0;
      79      108573 :   for (j=0;j<as->n->n;j++) {
      80      185699 :     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
      81       54288 :     if (as->x[0]->ops->mdot_local) PetscUseTypeMethod(as->x[j],mdot_local,n,bx,work);
      82           0 :     else PetscCall(VecMDot(as->x[j],n,bx,work));
      83      185699 :     for (i=0;i<n;i++) r[i] += work[i];
      84             :   }
      85             : 
      86             :   /* If def(__WITH_MPI__) and exists mdot_local */
      87             : #if defined(__WITH_MPI__)
      88       54285 :   if (as->x[0]->ops->mdot_local) {
      89             :     /* z[i] <- Allreduce(work[i]) */
      90      162855 :     PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
      91             :   }
      92             : #endif
      93             : 
      94       54285 :   PetscCall(PetscFree2(work0,bx));
      95       54285 :   PetscFunctionReturn(PETSC_SUCCESS);
      96             : }
      97             : 
      98           3 : static PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
      99             : {
     100           3 :   PetscScalar    sum = 0.0,work;
     101           3 :   PetscInt       i;
     102           3 :   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
     103             : 
     104           3 :   PetscFunctionBegin;
     105           3 :   SlepcValidVecComp(a,1);
     106           3 :   SlepcValidVecComp(b,2);
     107           3 :   if (as->x[0]->ops->tdot_local) {
     108           9 :     for (i=0,sum=0.0;i<as->n->n;i++) {
     109           6 :       PetscUseTypeMethod(as->x[i],tdot_local,bs->x[i],&work);
     110           6 :       sum += work;
     111             :     }
     112             : #if defined(__WITH_MPI__)
     113           3 :     work = sum;
     114           9 :     PetscCallMPI(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
     115             : #endif
     116             :   } else {
     117           0 :     for (i=0,sum=0.0;i<as->n->n;i++) {
     118           0 :       PetscCall(VecTDot(as->x[i],bs->x[i],&work));
     119           0 :       sum += work;
     120             :     }
     121             :   }
     122           3 :   *z = sum;
     123           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     124             : }
     125             : 
     126           3 : static PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
     127             : {
     128           3 :   PetscScalar    *work,*work0,*r;
     129           3 :   Vec_Comp       *as = (Vec_Comp*)a->data;
     130           3 :   Vec            *bx;
     131           3 :   PetscInt       i,j;
     132             : 
     133           3 :   PetscFunctionBegin;
     134           3 :   SlepcValidVecComp(a,1);
     135           9 :   for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
     136             : 
     137           3 :   if (as->n->n == 0) {
     138           0 :     *z = 0;
     139           0 :     PetscFunctionReturn(PETSC_SUCCESS);
     140             :   }
     141             : 
     142           3 :   PetscCall(PetscMalloc2(n,&work0,n,&bx));
     143             : 
     144             : #if defined(__WITH_MPI__)
     145           3 :   if (as->x[0]->ops->mtdot_local) {
     146           3 :     r = work0; work = z;
     147             :   } else
     148             : #endif
     149             :   {
     150           0 :     r = z; work = work0;
     151             :   }
     152             : 
     153             :   /* z[i] <- a.x' * b[i].x */
     154           9 :   for (i=0;i<n;i++) r[i] = 0.0;
     155           9 :   for (j=0;j<as->n->n;j++) {
     156          18 :     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
     157           6 :     if (as->x[0]->ops->mtdot_local) PetscUseTypeMethod(as->x[j],mtdot_local,n,bx,work);
     158           0 :     else PetscCall(VecMTDot(as->x[j],n,bx,work));
     159          18 :     for (i=0;i<n;i++) r[i] += work[i];
     160             :   }
     161             : 
     162             :   /* If def(__WITH_MPI__) and exists mtdot_local */
     163             : #if defined(__WITH_MPI__)
     164           3 :   if (as->x[0]->ops->mtdot_local) {
     165             :     /* z[i] <- Allreduce(work[i]) */
     166           9 :     PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
     167             :   }
     168             : #endif
     169             : 
     170           3 :   PetscCall(PetscFree2(work0,bx));
     171           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     172             : }
     173             : 
     174       57454 : static PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
     175             : {
     176       57454 :   PetscReal      work[3],s=0.0;
     177       57454 :   Vec_Comp       *as = (Vec_Comp*)a->data;
     178       57454 :   PetscInt       i;
     179             : 
     180       57454 :   PetscFunctionBegin;
     181       57454 :   SlepcValidVecComp(a,1);
     182             :   /* Initialize norm */
     183       57454 :   switch (t) {
     184          12 :     case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
     185       57439 :     case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
     186           3 :     case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
     187             :   }
     188      114932 :   for (i=0;i<as->n->n;i++) {
     189       57478 :     if (as->x[0]->ops->norm_local) PetscUseTypeMethod(as->x[i],norm_local,t,work);
     190           0 :     else PetscCall(VecNorm(as->x[i],t,work));
     191             :     /* norm+= work */
     192       57478 :     switch (t) {
     193           6 :       case NORM_1: *norm+= *work; break;
     194       57448 :       case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
     195           6 :       case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
     196          24 :       case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
     197             :     }
     198       57478 :   }
     199             : 
     200             :   /* If def(__WITH_MPI__) and exists norm_local */
     201             : #if defined(__WITH_MPI__)
     202       57454 :   if (as->x[0]->ops->norm_local) {
     203       57454 :     PetscReal work0[3];
     204             :     /* norm <- Allreduce(work) */
     205       57454 :     switch (t) {
     206           3 :     case NORM_1:
     207           3 :       work[0] = *norm;
     208           9 :       PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a)));
     209             :       break;
     210       57439 :     case NORM_2: case NORM_FROBENIUS:
     211       57439 :       work[0] = *norm; work[1] = s;
     212      172317 :       PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
     213       57439 :       *norm = GetNorm2(work0[0],work0[1]);
     214       57439 :       break;
     215           3 :     case NORM_1_AND_2:
     216           3 :       work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
     217           9 :       PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
     218           3 :       norm[0] = work0[0];
     219           3 :       norm[1] = GetNorm2(work0[1],work0[2]);
     220           3 :       break;
     221           9 :     case NORM_INFINITY:
     222           9 :       work[0] = *norm;
     223          27 :       PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a)));
     224             :       break;
     225             :     }
     226       57454 :   }
     227             : #else
     228             :   /* Norm correction */
     229           0 :   switch (t) {
     230           0 :     case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
     231           0 :     case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
     232           0 :     default: ;
     233             :   }
     234             : #endif
     235       57454 :   PetscFunctionReturn(PETSC_SUCCESS);
     236             : }
     237             : 
     238        1469 : PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
     239             : {
     240        1469 :   PetscScalar       dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
     241        1469 :   const PetscScalar *vx,*wx;
     242        1469 :   Vec_Comp          *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
     243        1469 :   PetscInt          i,n;
     244        1469 :   PetscBool         t0,t1;
     245             : #if defined(__WITH_MPI__)
     246         733 :   PetscScalar       work[4];
     247             : #endif
     248             : 
     249        1469 :   PetscFunctionBegin;
     250             :   /* Compute recursively the local part */
     251        1469 :   PetscCall(PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0));
     252        1469 :   PetscCall(PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1));
     253        1469 :   if (t0 && t1) {
     254         733 :     SlepcValidVecComp(v,1);
     255         733 :     SlepcValidVecComp(w,2);
     256        1469 :     for (i=0;i<vs->n->n;i++) {
     257         736 :       PetscCall(VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1));
     258         736 :       dp0 += dp1;
     259         736 :       nm0 += nm1;
     260             :     }
     261         736 :   } else if (!t0 && !t1) {
     262         736 :     PetscCall(VecGetLocalSize(v,&n));
     263         736 :     PetscCall(VecGetArrayRead(v,&vx));
     264         736 :     PetscCall(VecGetArrayRead(w,&wx));
     265       81052 :     for (i=0;i<n;i++) {
     266       80316 :       dp0 += vx[i]*PetscConj(wx[i]);
     267       80316 :       nm0 += wx[i]*PetscConj(wx[i]);
     268             :     }
     269         736 :     PetscCall(VecRestoreArrayRead(v,&vx));
     270         736 :     PetscCall(VecRestoreArrayRead(w,&wx));
     271           0 :   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");
     272             : 
     273             : #if defined(__WITH_MPI__)
     274             :     /* [dp, nm] <- Allreduce([dp0, nm0]) */
     275         733 :     work[0] = dp0; work[1] = nm0;
     276        2199 :     PetscCallMPI(MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v)));
     277         733 :     *dp = work[2]; *nm = work[3];
     278             : #else
     279         736 :     *dp = dp0; *nm = nm0;
     280             : #endif
     281        1469 :   PetscFunctionReturn(PETSC_SUCCESS);
     282             : }
     283             : 
     284             : #undef __SUF__
     285             : #undef __QUOTEME
     286             : #undef __SUF_C__