dtcsrmm.c

用于求解大型稀疏线性方程组Ax=b的数值计算库.

/* Call triangular mat-mult with upper triangular matrix: */

      transa = 1;
      dcsrmm( transa, m, n, k, alpha, descra, ua,
              uindx, upntrb, upntre, b, ldb,
              beta, c, ldc, work, lwork);

      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for upper(or lower) general matmult");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }
                                          /* Second, test symmetric matrices */
      printf("   Symmetric matrices:\n");

      descra[0] = 0;                      /* mat-mult with explicit matrix   */

      for (i=0;i!=m;i++)                  /* Initialize c                    */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call mat-mult with explicit symmtric matrix            */

      transa = 0;
      dcsrmm( transa, m, n, k, alpha, descra, a,
              indx, pntrb, pntre, b, ldb,
              beta, c, ldc, work, lwork);

      for (i=0;i!=n*m;i++)                /* copy result to d        */
        d[i] = c[i];

      descra[0] = 1;                      /* symmetry is implicit    */
      descra[1] = 1;                      /* lower triangular matrix */
      printf("      lower triangular\n");

      for (i=0;i!=m;i++)                  /* Initialize c            */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call symmetric mat-mult with lower triangular matrix     */

      transa = 0;
      dcsrmm( transa, m, n, k, alpha, descra, la,
              lindx, lpntrb, lpntre, b, ldb,
              beta, c, ldc, work, lwork);

                                          /* compare explicit to implicit */
      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for symmetric matmult (lower triangular)");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }

  
      descra[1] = 2;                      /* upper triangular matrix */
      printf("      upper triangular\n");

      for (i=0;i!=m;i++)                  /* Initialize c            */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call symmetric mat-mult with lower triangular matrix     */

      transa = 0;
      dcsrmm( transa, m, n, k, alpha, descra, ua,
              uindx, upntrb, upntre, b, ldb,
              beta, c, ldc, work, lwork);

                                          /* compare explicit to implicit */
      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for symmetric matmult (upper triangular)");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }
 
                                          /* Third, test skew-symmetric matrices */
      printf("   Skew-Symmetric matrices:\n");

      descra[0] = 0;                      /* mat-mult with explicit matrix   */

      for (i=0;i!=m;i++)                  /* Initialize c                    */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call mat-mult with explicit skew-symmetric matrix            */

      transa = 0;
      dcsrmm( transa, m, n, k, alpha, descra, ka,
              indx, pntrb, pntre, b, ldb,
              beta, c, ldc, work, lwork);

      for (i=0;i!=n*m;i++)                /* copy result to d        */
        d[i] = c[i];

      descra[0] = 4;                      /* symmetry is implicit    */
      descra[1] = 1;                      /* lower triangular matrix */
      printf("      lower triangular (no transp)\n");

      for (i=0;i!=m;i++)                  /* Initialize c            */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call skew-symmetric mat-mult with lower triangular matrix     */
      transa = 0;
      dcsrmm( transa, m, n, k, alpha, descra, la,
              lindx, lpntrb, lpntre, b, ldb,
              beta, c, ldc, work, lwork);

                                          /* compare explicit to implicit */
      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for skew-symmetric matmult (lower triangular)");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }

  
      descra[1] = 2;                      /* upper triangular matrix */
      printf("      upper triangular (transp)\n");

      for (i=0;i!=m;i++)                  /* Initialize c            */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call skew-symmetric mat-mult with upper triangular matrix     */
      transa = 1;                         /* Use transpose to get same    */
                                          /* matrix as lower triangular   */
      dcsrmm( transa, m, n, k, alpha, descra, ua,
              uindx, upntrb, upntre, b, ldb,
              beta, c, ldc, work, lwork);

                                          /* compare explicit to implicit */
      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for skew-symmetric matmult (upper triangular)");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }
                                          /* Now, work with transp of lower */
                                          /* and upper triangular matrix,   */
                                          /* results should be negation of  */
                                          /* explicit matrix multiply       */
                                          /* check by taking alpha = -alpha */

      descra[1] = 1;                      /* lower triangular matrix */
      printf("      lower triangular (transp)\n");

      for (i=0;i!=m;i++)                  /* Initialize c            */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call skew-symmetric mat-mult with lower triangular matrix     */
      transa = 1;
      dcsrmm( transa, m, n, k, -1.0*alpha, descra, la,
              lindx, lpntrb, lpntre, b, ldb,
              beta, c, ldc, work, lwork);

                                          /* compare explicit to implicit */
      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for skew-symmetric matmult (lower triangular)");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }

  
      descra[1] = 2;                      /* upper triangular matrix */
      printf("      upper triangular (no transp)\n");

      for (i=0;i!=m;i++)                  /* Initialize c            */
        for (j=0;j!=n;j++)
          c[j*m+i] = i+1;

/* Call skew-symmetric mat-mult with upper triangular matrix     */
      transa = 0;                         /* Use transpose to get same    */
                                          /* matrix as lower triangular   */
      dcsrmm( transa, m, n, k, -1.0*alpha, descra, ua,
              uindx, upntrb, upntre, b, ldb,
              beta, c, ldc, work, lwork);

                                          /* compare explicit to implicit */
      error = resid(n*m, d, c);
      if ( error >= tolerance ){
         errcount++;
         printf("Error for skew-symmetric matmult (upper triangular)");
         printf("n = %d.\n",n);
         printf("Residual: %10.6f \n",error);
         for (i=0;i!=n*m;i++) 
          printf("%6.2f  %6.2f\n",d[i], c[i]);
      }


} /* close loop on n */

if ( errcount > 0 ) 
  printf("%d errors in dtcsrmm run for alpha = %e, beta = %e\n",errcount,alpha, beta);

return errcount;
} /* end main */

double resid(int m, double *x1, double *x2) 
{

   double norm;
   int i;

   norm = 0.0;
 
   for (i=0;i<m;i++) norm += abs(x1[i] - x2[i]);

   if ( m == 0 ) {
     return norm;
   } else {
     return norm/m;
   }
 
}